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CONVERSION  TABLE 
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i 

Conversion  factors  for  U.S.  customary  to  metric  (SI)  units  of  measurement  ! 


To  Convert  Prom 

To 

Multiply 

angstrom 

neters  Im) 

1.000  000  X  £-10 

atmosphere  (normal) 

kilo  pascal  (kPa) 

1.013  25  XE+2 

bar 

kilo  pascal  (kPa) 

1.0C0  000  X  E+2 

barn 

meter2  (m2) 

1.000  000  XE-2S 

British  Thermal  unit  (thermochemical) 

Joule  (J) 

1.054  350  X  E+3 

calorie  (ttiermochemlcai) 

Joule  (J) 

4.184  000 

cal  (tbermochc-mlcnll/crn2 

mega  Joule  /  m*(MJ  /  m2) 

4.184  000  XE-2 

curie 

glga  becquerel  (GBq)* 

3.700  000  XE+l 

degree  (angle) 

radian  (rad) 

1.745  329  XE-2 

degree  Fahrenheit 

degree  kelvln  (K) 

tK»(t°f+  459.671/1.8 

electron  volt 

Joule  (J) 

1.602  19XE-19 

*** 

Joule  (J) 

1.000  000  X  E-7 

erg/ second 

watt  (W) 

1.000  000  XE-7 

foot 

meter  (ml 

3.048  000XE-1 

foot-pound-force 

joule  (J) 

1.355  818 

gallon  (U.S.  liquid) 

meter3  (m3) 

3.785  412  XE-3 

inch 

meter  (ml 

2.540  000  XE-2 

Jerk 

Joule  (J) 

1.000  000  XE+fi 

joule /kilogram  (J/Kg)  (radiation  dose 
absorbed! 

Gray  (Gy) 

1.000  000 

kllotons 

terajoules 

4.183 

kip  ( 1000  tbf) 

newton  (N) 

4.448  222  X  E+3 

kip /inch2  (ksi) 

kilo  pascal  (kPa) 

6.894  757  X  E+3 

ktap 

newton-second /m2  (N-s/m2) 

1.000  000  XE+2 

micron 

meter  (ml 

1.000  000  XE-6 

mil 

meter  (m) 

'  2.540  OOOXE-5 

mile  (international) 

meter  (ra| 

1.609  344  XE+3 

ounce 

kilogram  (kg) 

2.834  S52  X  E-2 

pound-force  (Ibf  avoirdupois) 

newton  (N) 

4.448  222 

pound-force  inch 

newton-meter  IN-m) 

1  129  848  X  E-l 

pound-force  /  inch 

newton/meter  (N/m) 

1.751  268  X  E+2 

pound-force  /  foot2 

kilo  pascal  (kPa) 

4.788  026  X  E-2 

pound-force /inch2  (p?l) 

kilo  pascal  (kPa) 

6.894  757 

pound-mass  (Ibm  avoirdupois) 

kilogram  (kg) 

4.535  924  X  E-l 

pound-mass-foot2  (moment  of  Inertia) 

kilogram-meter2  (kg-m2) 

4.214  011  XE-2 

pound-mass  /  foot3 

kilogram /meter3  (kg/m3) 

1.601  846  XE+l 

rad  (radiation  dose  absorbed) 

Gray  (Gy)** 

1.000  000  X  E-2 

roentgen 

coulomb /kilogram  (C/kg) 

2.579  760  X  E-4 

shake 

second  (s) 

1.000  000XE-8 

slug 

kilogram  (kg) 

1.459  390  XE+l 

torr  (mm  Hg.  0°C1 

kilo  pascal  (kPa) 

1.333  22  X  E-l 

The  becquerel  (Bq)  Is  the  SI  unit  of  radioactivity:  Bp  ■  1  eveni/s. 
••The  Gray  (Gy)  Is  the  SI  unit  of  absorbed  radiation. 
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SECTION  1 


INTRODUCTION 


1.1  BACKGROUND  AND  OBJECTIVE. 

In  1989  the  Defense  Nuclear  Agency  (DNA)  initiated  the  Underground  Technology 
Program  (UTF),  a  multi-year  investigation  into  the  vulnerability  of  underground  structures.  The 
program  includes  calculations,  material  modelling,  laboratory  testing,  and  field  testing  all  aimed 
at  improving  the  ability  to  predict  the  response  and  failure  of  underground  structures  subjected 
to  ground  shock  due  to  near-surface  explosions.  The  emphasis  is  on  deeply  buried  tunnels  with 
little  or  no  reinforcement  The  focal  point  of  the  program  is  a  field  test  to  be  conducted  at 
Ft  Knox,  wherein  a  buried  high-explosive  charge  will  be  used  to  load  various  deeply  buried 
structures  in  a  saturated  limestone  formation. 

During  the  course  of  the  UT?  at  least  five  different  organizations  have  been  engaged  in 
numerical  modelling  of  various  aspects  of  the  relevant  dynamic  processes.  Each  organization  has 
its  own  preconceptions  about  how  best  to  construct  these  models.  In  order  to  explore  the 
influence  of  the  particular  computational  approach  on  the  outcome  of  the  numerical  simulation, 
DNA  conducted  a  "benchmark  calculation"  exercise.  Several  idealized  problems  were  defined, 
and  each  participating  organization  was  asked  to  apply  its  computational  tools  to  generate  certain 
specified  outputs.  This  report  summarizes  the  results  of  that  activity. 

Large  rock  masses  inevitably  contain  discontinuities  which  are  referred  to  as  joints.  In 
many  cases  they  occur  in  one  or  several  parallel  sets  with  individual  joints  spaced  at  regular 
intervals.  The  numerical  treatment  of  the  combined  effect  of  motion  across  joints  and 
deformation  of  the  intervening  continuous  blocks  poses  a  formidable  challenge.  Most  of  the 
participants  use  two  general,  complementary  approaches:  explicit,  in  which  the  motions  across 
the  joints  and  the  deformations  within  the  blocks  are  represented  separately;  and  implicit,  wherein 
a  single,  effective  medium  is  defined  so  as  to  deform  on  the  average  like  an  assemblage  of  blocks 
and  joints1.  An  implicit  model,  once  defined,  can  be  (and  was  in  this  study)  used  just  like  any 
other  material  model  in  3  general  purpose  finite-element  or  finite  difference  code,  while  explicit 
models  require  special  treatment.  The  model  types  are  complementary  in  the  sense  that  when 
applied  to  tunnel  failure  problems,  the  implicit  niodels-which  are  more  computationally  efficient- 
suffice  in  regions  remote  from  the  tunnel,  while  the  explicit  ones-which  are  locally  mere 
accurate-must  be  used  near  the  tunnel  if  the  details  of  individual  block  motions  are  to  be 
captured.  As  expected,  the  approaches  to  joint  modelling  provided  the  most  stark  contrasts 
among  the  participants. 


*In  this  report  the  terms  explicit  and  implicit  will  always  refer  to  the  method  for  representing 
joints.  Tney  will  not  refer  to  the  time  integration  scheme,  except  in  the  following  sentence.  All 
participants  used  explicit  time  integration. 
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1.2  SCOPE 


Table  1-1  contains  a  summary  of  the  problems  posed  in  the  UTP  benchmark  calculation 
exercise.  More  precise  definitions  will  be  given  later,  when  each  problem  is  discussed  in  detail. 
Problems  1-IN  and  1-IM  exercise  the  intact  (unjointed)  and  implicitly  jointed  material 
respectively  over  the  same  simple  strain  path.  Problems  2-EX  and  2-IM  contrast  the  explicit  and 
implicit  treatment  of  a  jointed  sample  over  a  simple  strain  path.  Problem  2-S  is  the  first  one  to 
cause  shearing  displacement  across  the  joint.  It  was  solved  explicitly  by  most  but  implicitly  by 
one  participant.  Problem  3  is  the  first  one  where  dynamics  and  large-scale  spatial  variation  come 
into  play.  It  tests  compatibility  between  adjacent  regions  of  explicitly  and  implicitly  jointed 
material.  Problem  4  represents  the  way  a  "real"  tunnel  failure  problem  might  be  modelled,  and 
provides  a  platform  for  comparing  the  tuunei  deformations  computed  by  the  various  schemes. 


Table  1-1.  The  UTP  benchmark  problems. 


Problem 

Name 

Abbreviation 

Geometry 

Loading  f 

1-Intact 

1-IN 

Single  element,  no  joints 

Quasi-static  strain  path 
representing  spherically 
divergent  flow 

1-Implicit 

1-IM 

Single  element  with  two 
perpendicular  joints 

Same  strain  path  as  1-IN 

2-Explicit 

2-EX 

’’Stack  of  bricks,"  i.e., 
several  continuous 
horizontal  joints,  several 
staggered  vertical  joints 

Quasi-static  uniaxial  strain 
compression  and  unload  | 

I  2-Implicit 

2-IM 

Same  as  2-EX,  but 
modelled  implicitly 

Same  as  2-EX 

|  2-Shear 

2-S 

Two  triangular  blocks 
separated  by  3  joint 

Quasi-static  boundary 
displacements  consistent  with 
uniform  plane  strain,  designed 
to  cause  joint  slippage 

3 

3 

Wedge  section  of  annulus, 
modelled  implicitly, 
except  for  inner  region, 
modelled  explicitly 

Dynamic  radial  stress  pulse  on 
inner  boundary 

4 

2  ■  aa-  ..-ja  vi  Bsr 

4 

Same  as  3  but  lined  tunnel 
within  inner  region 

Same  as  3 

Table  1-2  lists  the  organizations  and  some  of  the  individuals  who  participated  in  the 
exercise.  It  also  contains  the  names  of  the  codes  used  to  perform  the  calculations.  The  only 


ones  which  are  generally  available  are  PRONTO  (through  Sandia  National  Laboratory)  and 
UDEC  (through  Itasca  for  a  fee);  the  others  are  proprietary. 
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Table  1-2.  Participants  in  the  UTP  benchmark  calculations!  exercise. 


Organization 

Abbreviation 

Point  of 
Contact 

Code  I 

California  Research  and  Technology 
Division,  The  Titan  Corporation 

CRT 

Y.  Marvin  Ito 

EXCAUBUR 

Itasca  Consulting  Group,  Inc. 

Itasca 

Loren  Lorig 

UDEC 

Lawrence  Livermore  National 
Laboratory 

LLNL 

Francois  Heuze 

DEBS 

RE/SPEC  Inc. 

RE/SPEC 

John  Osnes 

PRONTO 

Weidlinger  Associates 

WA 

Felix  W'ong 

FLEX  | 

1  Logicon  R  &  D  Associates 

RDA 

Don  Simons 

Exact  Solutions  I 

A  few  brief  remarks  about  the  various  approaches  may  be  of  interest  at  this  point.  All 
except  Itasca  and  LLNL  use  finite  elements  to  represent  the  rock  as  a  deformable  continuum; 

Itasca  uses  finite  difference  elements  while  LLNL  treats  the  rock  as  rigid  and  approximates  the 

rock’s  defonnability  by  modifying  the  stiffness  of  the  joints.  All  except  CRT  use  some  form  of 

a  slide  line  as  part  or  all  of  their  explicit  joint  model;  CRT  uses  a  special  finite  element.  CRT  ( 

has  the  most  sophisticated  implicit  model,  one  which  admits  arbitrary  constitutive  behavior  in 

both  the  intack  rock  and  joint,  and  by  enforcing  internal  compatibility  and  stress  equilibrium 

derives  a  super-element  representing  the  combined  deformation  due  to  both  joints  and  intact  rock. 

Wa  and  Itasca’s  implicit  models  are  isotropic  elastic-plastic,  while  CRT  and  RE/SPEC’s  are 

anisotropic.  Since  the  rock  in  this  study  contains  two  orthogonal  joint  sets,  the  large-scale  rock 

mass  behavior  will  be  orthotropic  (see  Section  4.3.3).  The  latter  two  implicit  modelling  * 

approaches  have  the  potential  to  represent  features  of  the  deformation  and  stress  fields  peculiar 

to  anisotropic  media. 

No  further  attempt  will  be  made  at  this  point  to  describe  the  details  of  these  codes’ 
treatments  of  the  benchmark  problems.  Instead,  each  organization  has  prepared  its  own  t 

description  which  is  included  in  this  report  as  an  appendix. 
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SECTION  2 


GENERAL  SPECIFICATIONS 


2.1  GEOMETRY. 

For  those  problems  involving  joints  (all  except  1-IN),  a  tunnel  (4),  or  complex 
computational  meshes  (3  and  4),  certain  geometric  parameters  are  standardized  as  shown  in 
Table  2-1.  Problems  3  and  4  are  in  plane  strain,  with  a  wedge-shaped  mesh  bounded  by  circular 
arcs  of  radii  450  and  550  m,  a  symmetry  line  through  the  tunnel  location,  and  a  ray  at  an  angle 
of  tan'10.1  =  5.71  (making  the  mesh  50  m  wide  at  the  tunnel  location).  An  inner,  rectangular 
region  extending  12.5  m  in  each  direction  from  the  tunnel  location  is  modelled  explicitly;  the 
rest,  implicitly. 


Table  2-1.  Standardized  geometric  parameters. 


Dimension 

Symbol 

Value  I 

Joint 

Joint  Spacing 

S 

1  m  | 

Joint  Thickness 

6 

5  mm  I 

Tunnel 

Range  of  Tunnel 

K 

500  m  | 

Tunnel  Diameter 

D 

5  m  I 

Liner  Thickness 

T 

50  mm 

|  Computational 

1  Mesh 

Discrete  Jointing  Zone 

L 

25m 

Mesh  Height 

H 

100  m 

Mesh  Divergence 

e 

tan*1  0.1 

2.2  MATERIAL  MODELS. 

The  intact  rock  is  treated  as  an  isotropic,  linear  elastic  and  perfectly  plastic  material,  with 
parameters  summarized  in  Table  2-2.  The  plastic  portion  has  an  associated  flow  rule  on  a  fixed 
(i.e.,  non-hardening)  two-invariant  stress  surface  fitted  to  a  Mohr-Coulomb  failure  condition  in 
triaxial  compression.  (This  is  sometimes  called  a  Mises-Schleicher  material.) 
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Table  2-2.  Material  properties. 


Property 

Symbol 

Value 

Intact  Mass  Density 

P 

2500  kg/m3 

Intact  Young’s  Modulus 

E 

30  GPa 

Intact  Poisson’s  Ratio 

v 

0.25 

Intact  Cohesion 

c 

4.5  MPa 

Intact  Friction  Angle 

<j> 

25° 

Intact  Dilation  Angle 

V 

25° 

Intact  Tensile  Strength 

T 

2  MPa 

Joint  Normal  Stiffness 

kN 

kN  -  kPa/m 

(b-iijf 

Joint  Shear  Stiffness 

kS 

1.25  GPa/m 

Joint  Cohesion 

CJ 

0 

Joint  Friction  Angle 

$ 

20°,  except  30°  in  problem  2-S 

Joint  Dilation  Angle 

W 

0 

Liner  Young’s  Modulus 

El 

200  GPa 

Liner  Poisson’s  Ratio 

VL 

0.30 

Liner  Yield  Strength 

400  MPa 

Liner  Mass  Density 

Pl 

7500  kg/m3 

To  be  more  specific,  the  failure  condition  in  terms  of  the  principal  stresses  a1?a2,a3 
(assumed  positive  in  compression)  can  be  written 

f(ova2,o3)  -  3J2  -( a+bp )2  -  0  ,  (2-1) 

/  2  2  2 

where  3 J2  -  °i  +a2 +a3  ~aia2  -a2a3  -o3o1  and  p  -  (01+a2+o3)/3  .  When  this  is  fitted  to 
a  Mohr-Coulomb  failure  criterion  in  a  state  of  triaxial  compression  (a2-a3<o1)  we  find  the 
relation  between  (a,b)  and  (<j>,c)  to  be 

6ccosd>  ,  6  sin* 

a  - - —  ,  b  - - -I-  ,  (2-2) 

3  -  sin4>  3  -  sin<j> 

For  the  parameters  given  in  Table  2-2  we  have  a=9A9  MPa,  fc=0.984. 


i 


> 


> 


► 


> 


5 


This  materia!  specification  is  adopted  for  its  computational  simplicity.  However,  it  is 
unphysical  in  that  it  leads  to  excessive  dilatancy  when  compared  with  the  response  of  real  rocks. 
This  needs  to  be  kept  in  mind  when  assessing  the  results  of  problem  4,  where — perhaps  contrary 
to  intuition,  although  not  to  certain  lab  and  field  observations — all  calculators  except  one  will  be 
seen  to  predict  that  the  tunnel  suffers  a  decrease  in  the  length  of  its  springline  (horizontal 
diameter).  The  one  exception — LLNL — did  not  model  the  rock  as  a  continuum  and  consequently 
could  not  possibly  represent  its  plastic  response  accurately. 

The  joint  behavior  is  also  assumed  to  be  clastic-plastic  with  a  constant  friction  angle. 
However,  the  normal  (opening/closing)  elastic  behavior  is  nonlinear,  and  the  inelastic  behavior 
is  non-dilatant.  The  liner  is  elastic-perfectiy  plastic  with  a  Mises  failure  condition  and  flow  rule. 
Joint  and  liner  property  values  are  also  listed  in  Table  2-2. 


rn 
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SECTION  3 

STATEMENTS  AND  RESULTS  FOR  PROBLEMS  WITH  ANALYTICAL  SOLUTIONS 


3.1  ANALYTIC  TREATMENT  OF  MATERIAL  RESPONSE. 

The  first  five  problems  (see  Table  1-1)  r»re  c;uas>static  ones  driven  by  boundary 
displacements  which  are  consistent  with  homogeneous  (uniform)  strain  throughout  the  region  of 
interest.  This  does  not  mean  that  the  actual  strains  will  be  uniform;  if  there  are  joints  then  actual 
strains  will  not  be  uniform.  But  in  fact  the  stresses  and  strains  in  ;he  intact  rock  and  joint 
material  will  separately  be  homogeneous  at  each  point  of  the  imposed  strain  path.  This  opens 
the  possibility  of  direct  analytical  solution  of  these  problems  for  comparison  with  numerical 
results. 


Another  way  of  viewing  the  situation  is  that  each  of  these  problems  vc-luces  to  nothing 
more  than  finding  the  response  of  a  single  implicit  element  around  a  specified  strain  path.  This 
is  strictly  a  material  response  question;  equations  of  motion  or  compatibility  among  elements  play 
no  role  whatsoever.  It  is  curious  to  note  that  only  one  of  the  participants  (CRT)  produced  single¬ 
element  solutions  to  all  of  the  first  five  problems. 

All  five  of  these  problems  can  easily  be  written  in  terms  of  principal 
stresses  ava2,a 3  and  strains  .  Also,  in  no  case  will  any  strain  exceed  three  percent, 

so  a  small-strain  analysis  is  sufficient.  As  usual,  the  strain  increments  are  decomposed  intc 
elastic  and  plastic  parts 

dz  j  -  dz^  +  dz ^  ,  (3-1) 

with  similar  forms  for  dz2  and  dz2  .  The  elastic  increments  follow  Hooke’s  law 

-v(rfo2vrfa3)j  ,  (3-2) 


with  similar  forms  for  the  other  two  components.  The  elastic  response  is  easily  invertible: 


dax  - 


E(1  -v) 

(l+vXl-2v) 


dz^  + 


-l~(dzf+dzf) 
1  -v 


(3-3) 


If  /(o1,o2,a3)<0,  then  the  entire  response  is  elastic  and  Eq.  (3-3)  fully  describes  the  intact  rock 
deformation.  Moreover,  since  the  coefficients  are  constant,  Eq.  (3-3)  can  be  integrated  by 
inspection  over  any  change  in  strain. 
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The  plastic  increments  follow  from  the  associated  flow  iaw  applied  to  the  yield 
function  (2-1): 


X[ao1-p(a2+a3)^] 


(3-4) 


where  X,  is  a  constant  (to  be  determined)  and  similar  formulas  obtain  for  the  other  components. 
In  (3-4)  the  constant  a,f5,y  are  given  by 


a  = 


P 


Y 


lab 

“ 


(3-5) 


with  ( a,b )  given  by  (2-2).  When  the  intact  rock  is  plastic,  total  strains  are  the  sum  of  (3-2)  and 
(3-4).  An  additional  condition  is  that  of  continued  yielding,  which  can  either  be  written  as 
Eq.  (2-1)  or 

df  =  J-Ldo,  +JjLda7  +_^_da,  =  0  , 

do1  1  da2  2  do3  3  1  ’ 

where  the  derivatives  are  given  by  formulas  like  the  square-bracketed  quantity  in  (3-4).  The 
three  principal  stress  increments  can  be  found  in  terms  of  principal  strain  increments  and  current 
stresses  by  solving  the  system  of  ten  equations  (3-1),  (3-3),  (3-4),  (3-6)  in  the  ten 
unknowns  dat ,  dtf\  dzf^,  X  .  We  will  defer  this  step  until  discussing  specific  problems. 

It  will  be  useful  to  represent  the  joint  normal  response  in  terms  of  stress.  According  to 
Table  2-2  the  joint  stiffness  is 

t  .  =  A 

N  ’  C3-7) 


> 


» 


> 


» 
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where  A=625  kPa-m  and  Uj  =  joint  closure.  This  can  be  rearranged  to  give 


do 


Aduj 
(6-“j)2  ’ 


» 


This  can  easily  be  integrated,  giving 


o 


Auj 

(6^7)6 


» 
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which  can  be  inverted  to  give 


> 


Uj  = 


b2o 
A  +ab 


The  stiffness  (3-7)  can  now  be  written  in  terms  of  stress: 


(3-8) 

» 


3.2  PROBLEM  1-IN:  SPHERICALLY  DIVERGENT  STRAIN  PATH  IN  INTACT  ROCK. 
3.2.1  Statement  of  the  Problem. 

An  intact  rock  sample  is  strained  uniformly  and  quasi-statically  along  a  strain  path  w:‘t' 
two  equal  principal  strains.  ‘This  represents  a  spherically  symmetric  deformation,  where  the  two 
hoop  strains  e0  are  equal  and  generally  different  from  the  spherical  radial  strain  tr.  The 
geometry  and  strain  path2  are  shown  in  Figure  3-1.  The  first  leg  is  nearly  uniaxial 
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» 
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Figure  3-1.  Geometry  and  loading  in  problems  1-IN  and  1-IM. 


2In  this  p'ot  as  well  as  throughout  this  report,  strains  and  stresses  are  taken  as  positive  in 
compression. 
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compression  and  represents  the  passage  of  a  shock  front.  During  the  second  !eg,  tensile  hoop 
strain  accumulates  due  to  outward  radial  displacement.  In  the  third  leg,  the  radial  strain  is  fully 
relieved  but  a  residual  tensile  hoop  strain  remains.  The  problem  is  to  find  the  stresses  at  each 
point  along  the  given  strain  path. 

3.2.2  Analytic  Solution. 


Let  o1~cr  and  o2=a3=oe  be  the  spherical  radial  and  hoop  stresses  respectively,  and 
similarly  for  strains.  From  (3-2),  the  elastic  response  is 


diC)  =  j(d°r~2vdoe)  ,  dzf  =  i[(l-v)do6-Wcg 


(3-9) 


with  the  inverse 


do  = 
r 


£[(1  -v)dz{‘]  +2vdz{‘]] 
(1  +v)(l  -2v) 


doQ  = 


E{dzf  +vdz^) 
(1  +v)(i  -2v) 


(3-10) 


In  the  elastic  regime  these  equations  govern  the  full  response. 
From  (3-5)  the  plastic  strain  increments  must  obey 


dzf>  -  2X|l  -J 


f1  »  dtf  =  -xf1+^-](ar-°e)  •  (3-11) 

\  5 )  \  5  ) 


The  differential  form  of  the  continuing  yield  condition  (3-6)  reduces  to 


(3-12) 


Equations  (3-9),  (3-11),  and  (3-12)  combine  with  the  two  independent  forms  of  (3-1)  to  provide 
seven  linear  equations  in  the  unknowns  (dz^,  dt.^,  dz^’ ,  dzg\  dar,  dae,X)  .  They  can  be 
solved  for  the  stress  increments  during  plastic  loading,  giving 


da,  = 


£(3+2b)[(3  +2b)dzr+2(3  -b)dzQ] 
27(1  -2v)  +6b2(l  +v) 


E(3-b)[(3+2b)dzr+2(3-b)dzQ] 

“ae  = - 


(3-13) 


27(1  -2v)  +6bz(l  +v) 


Note  that  with  the  given  yield  function  and  flow  rule,  not  only  the  elastic  formula  (3-10) 
but  the  elastic-plastic  one  (3-13)  as  well  are  linear  with  constant  coefficients.  Thus  as  long 
as  dzQ/dzr  is  constant,  the  finite  stress  increments  obey  the  same  equations  as  the  infinitesimals. 

To  compute  the  evolution  of  stresses  along  the  strain  path  of  Figure  3-1,  we  must  proceed 
in  steps,  beginning  with  elastic  response.  On  the  initial  leg  OA  we  have  dzQ/dzr  =  -1/30  . 


•  •  • 
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The  stresses  at  any  point  along  the  leg  are  those  at  the  beginning — 0  in  this  case — added  to  the 
increments.  Combining  this  strain  increment  ratio  with  (3-10)  we  find 


da. 


e  _  30v-l 


da. 


30-32v 


=  0.2954 


(3-14) 


This  holds  until  yielding,  and  in  this  case  gives  the  ratio  of  stresses  themselves  since  they  both 
started  at  0.  For  this  axisymmetric  stress  state  the  yield  condition  (2-1)  can  be  recast 


/  =  |or-oe| 


+|(2or+oe) 


=  0  . 


(3-15) 


Equation  (3-14)  and  (3-15)  solved  simultaneously  give  the  stresses  at  first  yield,  labeled  point 
A1  in  Table  3-1.  The  strains  at  this  point  follow  from  (3-9). 


Table  3-1.  Analytic  solution  to  problem  1-IN. 


Point 

Time 

Radial 

Stress 

°r 

(MPa) 

Hoop 

Stress 

°0 

(MPa) 

Pressure 

P 

(MPa) 

Stress 

Difference 

0r-°8 

(MPa) 

Radial 

Strain 

Er 

(percent) 

Hoop 

Strain 

Ee 

(percent) 

Volume 

Strain 

(percent) 

0 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

A1 

0.49 

51.9 

15.3 

27.5 

36.6 

0.148 

-0.005 

0.138 

H 

1. 

104.8 

26.8 

59.5 

68.0 

0.3 

-0.01 

0.28 

B 

2. 

62.5 

19.6 

33.9 

42.9 

0.23 

-0.07 

0.09 

Bl 

2.59 

9.2 

-2. 

1.7 

11.2 

0.095 

-0.088 

-0.08 

B2 

2.59 

10.2 

0. 

mm 

10.2 

0.095 

-0.088 

-0.08 

B3 

2.76 

-2. 

0. 

-0.7 

-2. 

0.054 

-0.093 

-0.13 

B4 

2.76 

0. 

0. 

0. 

0. 

0.054 

-0.093 

-0.13 

C 

3. 

0. 

0. 

0. 

0. 

0. 

-0.1 

-0.2 

> 
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Once  yielding  occurs,  equations  (3-13)  obtain.  These  two  equations  together  with  the 
strain-path  slope  z^Jtr  =  -1/30  give  the  remaining  response  on  leg  OA  as  a  function  of  I 

(say)  dEr  and  the  final  value  er  -  0.3  defines  point  A. 


> 


» 


•  • 


When  the  strain  path  direction  changes  at  point  A,  we  can  tell  if  the  next  increment  is 
elastic  or  plastic  by  supposing  that  it  is  elastic  ("taking  a  trial  elastic  increment"),  in  which  case 


from  (3-10)  with  <2eGAier=6/7  we  find  daQ/dar  =31/33  .  With  o^>aQ  and  fc=0.984  the 
increment  in  the  yield  function  (3-15)  would  be 

df  =  dar  1  -  —  -  —  fl  + 1)  =  -0.904 do,  .  (3-16) 

r  3  3j 

With  dar< 0  the  trial  increment  would  cause  /to  go  positive,  which  is  prohibited;  therefore  leg 
AB  must  start  with  an  elastic-plastic  increment.  Starting  from  the  known  state  at  A,  as  before 
we  use  (3-13)  and  the  strain-path  slope  de0/der= 6/7  to  find  the  changes  in  stresses  as  a  function 
of  (say)  dtr .  Provided  the  tensile  failure  criterion  is  not  met,  the  response  will  continue  in  the 
elastic-plastic  mode  until  the  next  change  in  strain-path  direction. 

The  tensile  limit  is  given  as  2  MPa,  but  no  other  details  of  the  tensile  failure  law  are 
specified.  In  this  analysis  we  assume  this  limit  applies  separately  to  each  principal  stress,  and 
once  it  is  exceeded,  that  stress  is  set  to  zero. 

Tensile  failure  therefore  does  not  occur  before  reaching  point  B,  so  we  can  fill  in  the  state 
at  B  according  to  elastic-plastic  response. 

The  beginning  of  leg  BC  is  treated  the  same  as  was  AB,  and  is  likewise  found  to  be 
elastic-plastic.  As  before,  the  strain  and  stress  increment  ratios  give  the  direction  of  the  stress 
path,  and  the  starting  point  is  known  from  the  previous  leg.  The  hoop  stress  is  the  first  to  reach 
the  tensile  limit.  Here  the  hoop  stress  changes  abruptly  from  -2MPa  (point  Bl)  to  0  (point  B2), 
while  the  radial  strain  remains  fixed.  The  macroscopic  hoop  strain  also  stays  fixed,  but  after 
cracking,  it  will  contain  a  portion  due  to  the  crack,  while  the  intact  material  will  contract  in  the 
hoop  directions.  Assuming  the  stress  adjustments  during  crack  formation  are  elastic,  from  (3-9)  j 
with  dzr  =  0  and  daQ  =  +2  MPa  we  find  dar  - 1  MPa .  This  completes  the  state  at  B2.  (A 
yield  function  check  verifies  that  the  stress  adjustment  was  indeed  elastic.) 

The  next  leg  is  (as  a  trial)  taken  to  be  elastic,  with  o0« 0 .  Equation  (3-9)  then 
giving  dor  =Edtr.  Radial  tensile  failure  occurs  as  or  reaches  -2  MPa  at  point  B3.  Now  the 
radial  stress  also  drops  to  0  (point  B4)  with  no  change  in  macroscopic  strain,  and  remains  there 
until  the  end  of  the  strain  path  at  point  C. 

3.2.3  Numerical  Solutions. 

The  compressibility  and  stress  paths  from  the  various  numerical  solutions  to  problem  1-IN 
are  shown  in  Figures  3-2  and  3-3  respectively.  The  most  obvious  feature  of  these  comparisons 
is  the  inadequacy  of  the  LLNL  approach  for  representing  stresses  in  intact  rock.  It  has  the  wrong 
stiffness,  volumetric  hysteresis,  and  stress  path  over  the  full  duration  of  the  problem.  Among  the 
other  curve?,  up  to  the  point  of  tensile  failure,  all  agree  well  except  for  slightly  low  unloading 
pressures  in  the  CRT  model.  The  final  leg  of  the  compressibility  curve  suggests  that  all  the 
numerical  models  except  RE/SFEC’s  used  a  criterion  based  on  pressure  rather  than  principal 
stress,  and  that  after  tensile  failure  some  held  the  pressure  at  the  cutoff  value  rather  than  fixing 
it  at  0  .  Differences  in  intact  rock  tensile  failure  models  are  probably  of  negligible  consequence 
to  tunnel  failure  in  jointed  media. 


3.3  PROBLEM  1-IM:  SPHERICALLY  DIVERGENT  STRAIN  PATH  IN  IMPLICITLY 
JOINTED  ROCK. 

3.3.1  Statement  of  the  Problem. 

The  only  difference  between  this  problem  and  the  previous  one,  problem  1-IN,  is  the 
material,  which  now  contains  two  perpendicular  joint  sets  modelled  implicitly.  One  joint  set  is 
assumed  to  be  normal  to  the  r-direction  (horizontal  in  Figure  3-1).  Joint  properties  are  given  in 
Tables  2-1  and  2-2. 


3.3.2  Analytic  Solution. 

The  introduction  of  implicit  jointing  profoundly  affects  the  solution  in  several  ways:  by 
compounding  the  rock  and  joint  response  in  each  individual  direction,  by  introducing  nonlinearity 
on  account  of  the  joint’s  nonlinear  normal  stiffness,  and  by  upsetting  the  axisymmetry  on  account 
of  the  effective  anisotropy  induced  by  the  joints.  The  principal  stresses  still  align  with  coordinate 
directions  and  joints,  so  there  are  still  no  shear  stresses  or  shear  displacements  across  joints,  and 
it  still  suffices  to  consider  only  the  normal  stresses,  strains,  and  displacements. 


For  this  discussion,  let  the  z-axis  align  with  the  r-direction  and  the  normal  to  one  joint 
set,  and  the  x-axis  with  the  other  joint  set.  Consider  the  displacement  increments  for  a  single 
joint  and  intervening  intact  rock  (i.e.,  a  1-m  cube): 


dux  =  dux  +  dux  -  wjdtx  + 

duy  =  duy  =  wydty  ’ 
duz  =  duz  +  duz  =  + 


(3-17) 


where  superscripts  /,  J  refer  to  intact  material  and  joints  respectively,  and  wx  s=  wy  =  w.  =  1  m. 
The  joint  thickness  has  been  neglected  compared  with  the  block  width  in  relating  intact  rock 
strains  to  the  corresponding  displacements.  Recall  that  the  joint  stiffnesses  kx  ,  kz  depend 
respectively  on  the  stresses  ox  ,  oz  according  to  (3-8).  Combining  (3-17)  with  the  elastic 
constitutive  law  (3-2)  gives 


(3-18) 


\L  E  ' 

-v 

■41 

1 

Art 

Edux 

""  j 

V 

dox 

"'Ai 

X 

\  / 

Edu 

: 

1 

X 

-v 

do„ 

y\ 

y 

L  E  \ 

y 

symmetric 

— 

1+ - 7 

Eduz 

■ 

"*1 

d°z 

wz 

This  system  of  three  linear  equations  must  be  solved  simultaneously  to  get  the  elastic  stress 
increments  corresponding  to  given  displacement  increments  dux ,  duy ,  duz  .  The  system  as  a 
whole  is  nonlinear,  because  the  joint  stiffnesses,  and  consequently  the  coefficients  in  (3-18), 
depend  on  the  current  stress  level. 


Three  of  the  equations  governing  stress  increments  in  the  plastic  regime  are  found  by 
combining  the  elastic  strain  increments  (3-2)  with  the  plastic  ones  (3-4)  and  substituting  the  sum 
into  (3-17).  A  new  unknown  X  has  been  introduced,  and  the  fourth  and  last  required  equation 
is  the  condition  of  continuing  plastic  yielding  (3-6).  The  partial  derivatives  take  the  form  of  the 
square-bracketed  quantity  in  (3-4).  The  resulting  system  of  four  equations  in  four  unknowns  is 
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(3-19) 


where  in  the  axisymmctric  problem  at  hand  we  have  dux/wx  =  auy/wy  =  dtQ ,  duz/wz  -  dtf  . 
Note  that  the  matrix  in  (3-19)  is  symmetric.  This  is  a  consequence  of  the  associated  plastic  flow 
law. 


Either  system  (3-18)  or  (3-19)  could  be  inverted  in  closed  form  if  desired.  However,  the 
stress  dependence  of  the  resulting  expressions  for  the  stress  increments  would  render  them  very 
difficult  if  not  impossible  to  integrate  in  closed  form.  Therefore  it  was  decided  to  perform  both 
the  inversion  and  the  integration  numerically.  The  procedure  for  each  time  step  is  as  follows: 


•  •  • 


•  •  • 
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(1)  Solve  the  elastic  equations  (3-18),  using  current  stresses,  for  trial  stress  increments 
due  to  current  displacement  increments.  Note  that  overall  strain  increments  in  the 
problem  specification  (Figure  3-1),  when  applied  to  the  1-m  cube,  are  numerically 
equal  to  the  displacements  in  meters  appearing  on  the  right-hand  sides  of 
equations  (3-18)  or  (3-19).  A  crude  predictor-corrector  is  effected  by  updating  the 
stresses  and  then  re-solving  the  equations. 

(2)  Check  whether  the  yield  criterion  is  violated  by  the  current  stresses  augmented  by 
the  trial  increments.  If  not,  the  current  increment  is  elastic  and  the  computed  trial 
values  are  retained  as  actual  values. 

(3)  If  the  current  state  is  elastic  but  the  trial  values  violate  the  yield  criterion,  linearly 
interpolate  back  to  the  place  where  f=  0.  reset  the  time  axis  to  this  point,  and  move 
forward  from  there  with  a  plastic  increment  based  on  equations  (3-19). 

(4)  If  the  current  state  is  plastic,  first  conduct  an  elastic  trial  as  above.  If  the  trial  state 
is  elastic,  update  immediately  (no  interpolation  is  needed);  if  not,  compute  a 
plastic  increment  from  (3-19)  and  update  immediately. 

It  was  not  deemed  of  sufficient  interest  to  pursue  this  computation  beyond  tensile  failure, 
so  it  was  simply  terminated  when  any  principal  stress  reached  the  tensile  limit.  The  results  are 
given  with  the  numerical  solutions  in  the  following  section. 

3.3.3  Numerical  Solutions. 

Only  CRT  and  RE/SPEC  submitted  numerical  solutions  to  this  problem.  Because  it  is  not 
axisymmetric,  the  stress  difference  oz~ax  is  not  equal  to  (37^) 1/2  ,  so  three  plots  of  results 
are  given.  Figures  (3-4),  (3-5),  and  (3-6)  indicate  that  prior  to  tensile  failure,  both  numerical 
results  are  quite  close  to  the  analytical  one.  After  that,  although  the  analytic  solution  has  not 
been  carried  this  far,  CRT’s  compressibility  curve  cannot  be  correct,  since  it  does  not  terminate 
at  a  volumetric  strain  of  -0.2  percent.  However,  as  noted,  this  portion  of  the  curve  is  not  of 
particular  interest  in  this  problem. 

3.4  PROBLEMS  2-EX  AND  2-IM:  UNIAXIAL  STRAIN  OF  A  "STACK  OF  BRICKS", 
3.4.1  Statement  of  the  Problem. 

A  4-by-5-by-l-m  volume  is  filled  with  1-m  cubes  of  rock  separated  by  joints  as  shown 
in  Figure  3-7.  There  are  no  joints  parallel  to  the  plane  of  the  figure.  Five  sides  have  roller 
boundaries.  The  top  is  displaced  vertically  and  quasi-statically  as  indicated  in  the  figure,  while 
shear  traction  is  held  at  zero.  Half-thickness  joints  separate  full  blocks  from  the  three  roller 
boundaries  normal  to  the  plane  of  the  figure. 


•  •  • 


•  •  •  • 
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Note:  Plane  strain  in  third  dimension 
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Figure  3-7.  Geometry  and  loading  in  problems  2-EX  and  2-IM. 


3.4.2  Analytical  Solution. 


No  distinction  between  implicit  and  explicit  rock  models  needs  to  be  made  for  the  purpose 
of  deriving  the  analytic  solution  to  this  problem.  The  stated  conditions  lead  to  stresses  and 
strains  identical  to  those  in  an  infinite  region  with  the  same  joint  spacing  and  homogeneous 
macroscopic  strain.  The  principal  directions  of  the  geometry  will  coincide  with  principal 
directions  of  stress  and  strain.  Equations  (3-18)  and  (3-19)  from  the  previous  section  govern  a 
single,  representative  1-m  cubic  element  under  these  conditions.  (The  staggering  in  the  jointing 
only  affects  deformations  at  a  much  higher  level  of  accuracy  than  considered  here.)  Therefore, 
the  analytic  solution  to  this  problem  follows  exactly  the  same  procedure  as  problem  1-IM.  The 
only  difference  is  the  right-hand  side  of  (3-18)  or  (3-19),  where  now  dux  =  du  =  0  ,  and  (for 
a  single  1-m  element)  uz  will  increase  quasistaticaily  up  to  0.3  cm  and  then  decrease  to  zero. 
The  results  of  the  step-by-step  numerical  integration  are  shown  in  Figures  3-8  and  3-9. 

3.4.3  Numerical  Solutions. 

Figures  3-8  and  3-9  also  show  the  compressibility  and  stress  path  according  to  the  various 
explicit  numerical  joint  models.  In  each  plot  two  curves  stand  out.  The  LLNL  model  again  has 
incorrect  stiffness,  hysteresis,  and  deviator  stresses.  It  does  appear  to  agree  with  the  initial  slope 
of  the  compression  loading  curve,  but  not  to  stiffen  on  loading  as  it  should  due  to  the  joint’s 
nonlinear  elastic  compressive  behavior. 

The  WA  model  has  incorrect  stiffness  and  excessive  numerical  chatter.  The  reason  for 
the  former  is  that  in  the  joint  treatment  used  in  this  problem,  joint  stiffness  was  fixed  solely  by 
numerical  considerations  related  to  the  detection  of  and  compensation  for  node  penetration  across 
joints.  There  was  no  other  mechanism  to  include  a  specified  normal  joint  stiffness,  so  it  was 
ignored.  The  required  numerical  stiffness  appears  to  have  been  much  greater  than  that  specified 
for  the  joint,  leading  to  stiffer  overall  response.  That  the  WA  stiffness  is  roughly  twice  the 
correct  value  is  consistent  with  the  fact  that  at  least  according  to  specifications  in  the  elastic 
regime,  the  contributions  of  the  intact  rock  and  joints  to  the  overall  stiffness  are  about  equal. 
Subsequent  to  this  calculation,  the  joint  model  was  modified  to  admit  a  specified  normal  elastic 
stiffness  in  place  of  the  former,  numerically  driven  value. 

The  excessive  chatter  in  both  compressibility  and  stress  path  in  the  WA  calculation  is  also 
caused  by  the  joint  model.  WA  explicitly  modelled  every  individual  block  and  interface  shown 
in  Figure  3-7.  Subsequent  to  this  problem,  WA  modified  the  treatment  of  joint  forces  at  the 
comers  of  blocks  and  believes  that  these  modifications  are  responsible  for  the  improvement  in 
the  model’s  behavior  in  subsequent  problems.  (Appendix  E  contains  a  brief  description  of  the 
modelling  approach  and  the  modifications.) 

Numerical  results  for  problem  2-IM  are  shown  in  Figures  3-10  and  3-11  together  with  the 
analytical  solution.  Again,  only  CRT  and  RE/SPEC  submitted  solutions  based  on  implicit  joint 
models,  and  both  agree  closely  with  the  analytical  results. 


3.5  PROBLEM  2-S:  SHEARING  OF  TWO  TRIANGULAR  BLOCKS  SEPARATED  BY  A 
SINGLE  JOINT. 

3.5.1  Statement  of  the  Problem. 

Since  none  of  the  foregoing  problems  exercised  the  joint  in  shear,  this  problem  was 
devised  to  do  just  that,  and  to  do  so  while  maintaining  a  homogeneous  strain  state  in  the  intact 
rock.  Two  blocks  cf  intact  rock  are  separated  by  a  diagonal  joint  as  shown  in  Figure  3-12. 
Plane  strain  applies  out  of  the  plane  of  the  figure.  Normal  displacements  are  imposed  on  the 
upper  and  right-hand  faces,  while  maintaining  zero  shear  traction.  The  remaining  boundaries  are 
rollers.  As  a  function  of  time  in  arbitrary  units,  the  displacements  are  given  in  Table  3-2,  with 
linear  variations  between  tabulated  values.  The  displacements  follow  the  path  indicated  in  the 
figure.  Tiie  objective  is  to  find  the  stress  histories  in  the  block.  The  problem  is  to  be  solved 
with  an  explicit  joint  model,  but  an  implicit  model  may  be  tried  in  addition. 


Figure  3-12.  Geometry  and  loading  in  problem  2-S. 


Table  3-2.  Imposed  boundary  displacements  in  problem  2-S. 


Time 

(Arbitrary 

units) 

“l 

(m) 

“2 

(m) 

0 

0 

0 

1 

0.0001 

0.0001 

2 

0.0187 

-0.0153 

3 

0.0175 

-0.0165 

4 

0.0005 

0.0005 

5 

0 

0 

3.5.2  Analytical  Solution. 

3.5.2.1  Outline  of  the  solution  procedure.  The  problem  will  be  treated  by  an 
incrementally  linear  scheme.  At  the  beginning  of  each  time  step  the  stresses,  boundary 
displacements,  total  joint  slip,  and  pending  displacement  increments  are  known,  and  the  stress 
increments  and  total  joint  slip  increment  are  to  be  found.  Once  the  increments  are  found,  the 
corresponding  quantities  are  updated  and  the  process  repeated. 

For  quasi-static  deformations,  this  problem  can  be  put  into  the  same  form  as  the  more 
conventional  one  of  elastic-plastic  deformation  with  two  intersecting,  non-hardening  failure 
surfaces.  One  surface  is  the  standard  one  for  the  intact  material;  the  other  is  the  locus  in  stress- 
space  of  points  satisfying  the  conditions  for  Mohr-Coulomb  frictional  slip.  Plastic  flow  in  intact 
rock  is  associated,  while  the  effective  plastic  strain  due  to  joint  slip  is  not  associated  with  the 
joint  failure  function.  Table  3-3  lists  the  four  potential  deformation  modes,  and  assigns  a  name 
and  number  to  each. 

Table  3-3.  Potential  deformation  modes  in  problem  2-S. 


Mode  number 

Mode  name 

Intact  response 

Joint  response 

1 

EE 

Elastic 

Elastic  (not  slipping) 

2 

EP 

Elastic 

Plastic  (slipping) 

3 

PE 

Plastic 

Elastic 

4 

PP 

Plastic 

Plastic 

•  •  • 


•  •  • 


•  • 
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For  each  potential  mode,  a  corresponding  set  of  stress  increments  can  always  be  obtained. 
The  formulas  will  be  given  below.  If  the  problem  is  well-posed,  then  one  and  only  one  of  these 
four  sets  will  be  admissible,  because  together  they  comprise  all  of  the  possibilities.  The  criteria 
for  admissibility  are  that  the  stresses  at  the  end  of  the  increment  be  inside  or  on  the  appropriate 
failure  surface(s),  and  that  the  plastic  strain  increments  for  the  active  plastic  mode(s)  be  in  the 
correct  direction.  The  former  requirement  is  equivalent  to  the  negativity  or  vanishing  of  the 
failure  function  / 1  or  fJ  .  For  the  intact  rock,  the  latter  requirement  is  that  the  proportionality 
constant  X 1  relating  plastic  strain  increments  to  stress  derivatives  of  the  flow  potential,  be 
positive.  For  the  joint,  the  condition  is  that  the  increment  of  slip  must  be  in  the  same  direction 
as  the  total  shear  stress  across  the  joint.  While  no  analytical  proof  of  uniqueness  is  available, 
it  will  be  found  numerically  for  the  prescribed  boundary  displacement  history  that  at  every  time 
step  there  is  one  and  only  one  admissible  solution  among  the  four  alternatives. 

The  computational  procedure  for  each  time  step  will  begin  with  the  evaluation  of  a 
number  of  "trial"  increments,  the  number  and  type  depending  on  the  location  of  the  stress  point 
at  the  start  of  the  increment.3  Specifically,  only  increments  of  types  which  cannot  be  ruled  out 
a  priori,  based  on  the  starting  point,  will  be  considered.  For  example,  if  the  starting  point  is  EE, 
i.e.,  both  failure  functions  are  negative  and  the  stress  point  is  inside  both  failure  surfaces,  then 
only  an  EE  increment  need  be  considered.  If  starting  with  an  EP  stress  state,  where  the  stress 
point  is  on  the  joint  failure  surface  but  inside  the  one  for  intact  material,  the  increment  must  be 
either  EP  or  EE.  If  a  trial  increment  causes  the  stress  point  to  cross  a  failure  surface,  then  linear 
interpolation  is  used  to  shorten  the  current  time  step  such  that  the  stress  point  will  be  right  on 
the  failure  surface  at  the  end  of  the  time  step.  This  ensures  that  only  a  single  deformation  mode 
occurs  during  the  entire  time  step. 

Table  3-4  is  a  road  map  of  the  entire  computational  scheme.  The  second  column  lists  the 
types  of  trial  increments  that  cannot  be  ruled  out  a  priori  for  each  starting  stress  point  location 
(column  1)  and  are  therefore  computed  and  checked  for  admissibility.  Conditions  for 
admissibility  are  listed  in  the  third  column  in  terms  of  failure  functions  and  plastic  flow 
proportionality  constants.  (Although  we  will  not  explicitly  write  the  joint  slip  increment  in  the 
form  of  a  constant  multiplying  the  gradient  of  a  flow  potential,  it  could  be  done.  In  the  table, 
"  X7*  >  0  "  is  simply  shorthand  for  "the  slip  increment  is  in  the  same  direction  as  the  total  shear 
stress  across  the  joint.")  The  superscript  refers  to  a  trial  value,  so  for 

example  f1*  =f,(p+do*)  where  o  is  the  stress  at  the  start  of  the  time  step  and  da *  is  the 
trial  stress  increment.  The  fourth  column  indicates  the  final  resolution  for  the  time  step,  based 
on  which  pair  of  admissibility  conditions  holds.  In  some  cases  the  trial  increment  itself  becomes 
the  actual  stress  increment.  In  others,  where  the  trial  increment  caused  the  stress  point  to  cross 
a  failure  surface,  an  interpolation  based  on  failure  function  values  is  applied  to  scale  back  both 
the  stress  increment  and  the  time  step  such  that  the  end  point  is  just  on  the  failure  surface(s). 


•^These  "trial"  increments  should  not  be  confused  with  the  trial  elastic  increment  used  in  most 
explicit  dynamic  codes. 
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Table  3-4.  Road  map  for  computational  scheme  used  in  problem  2-S. 


Stress  state  at 
start  of  time 
step 

Trial  increment 
types 

Admissible 
outcomes  (see 
note  A) 

Action  if 
condition  in 
column  3  holds 

Stress  state 
at  end  of 
time  step 

EE 

EE 

f!*<0JJ9<0 

da=da' 

EE 

/7*<0,  fJ*>0 

do=tfda * 

EP 

/7*>0,/J*<0 

da^tfda* 

PE 

/7*>0,  fJ*>0 

(See  note  B) 

EP  or  PE 

EP 

EE 

/7*<0,  fJ*<0 

da=da 

EE 

/7*>0,  fJ-<0 

do=Qtda* 

PE 

EP 

/7*< 0,  X/*>0 

da=da * 

EP 

/7*> o,  x/*>o 

da-^fda* 

PP 

PE 

EE 

/7*<0,  fJ*<0 

do=da * 

EE 

/7*<0,  fJ*>0 

da=tfda* 

EP 

PE 

)!*>0,fJ*<0 

da=da * 

PE 

X7*>0,/7*>0 

da-^da* 

PP 

PP 

EE 

/7*<0,  fJ\ 0 

do=do* 

EE 

EP 

/7*< 0,  kJ*>0 

da=da * 

EP 

PE 

X7*>0,  fJ*<0 

da=da ’ 

PE 

PP 

X7*>0,  */*>o 

da=da * 

PP 

Note  A:  For  uniqueness,  only  one  of  the  four  conditions  will  hold. 

This  condition,  corresponding  to  an  EE  starting  point  very  near  the 
intersection  and  a  strain  increment  outside  both  surfaces,  never  occurred  but 
could  be  handled  if  necessary. 


Note  B: 


I 


For  example,  the  constant  Q1 ,  which  must  fall  in  the  range  (0,1),  is  defined  by 


of  .  -f'W _  . 

f‘(p  .da*) - f'(a ) 

Finally,  the  last  column  shows  the  location  of  the  stress  point  at  the  end  of  the  time  step. 


3.5.22  Analytical  preliminaries.  As  before,  we  decompose  the  displacements  into  parts 
separately  due  to  intact  material  and  the  joint: 


“x  -  ux  +ux  *  *** 


uz  =  Uz  +  uz  =  M*z+«z 


u  =  u  =  0 
y  y 


(3-20) 


where  the  last  relation  expresses  the  plane  strain  condition  in  the 
absence  of  any  joints  normal  to  the  y-axis.  As  indicated  by  the 
vector  diagram  in  Figure  3-13,  the  parts  of  the  total  displacement 
due  to  motion  across  the  joint  can  be  further  decomposed  into 
components  normal  and  tangential  to  the  joint  These  are 
respectively  denoted  (positive  in  compression) 
and  u £  (positive  when  the  upper  block  moves  down  and  to  the 
right),  and  are  given  by 


J  J 
J  _  “n+«s 
uz 


& 


& 


(3-21)  figure  3-13.  Decomposition 
of  joint  displacements 


The  symmetries  of  this  problem  guarantee  that  the  only  nonvanishing  stresses  in  the  x-y-z 
coordinate  system  are  the  normal  stresses  in  the  coordinate  directions.  By  a  45-degree  cw 
coordinate  rotation,  these  give  rise  to  the  following  normal  and  shear  stresses  across  the  joint: 


» 


» 


> 


J  °z+°x  -  j  °z~ax  A  O 

°N  =  — 2  *  ’  5  "  — 2  "  T 


(3-22) 


In  view  of  the  form  of  (3-22)  and  other  factors  that  lead  to  later  simplification,  we  will 
regard  the  fundamental  stress  unknowns  to  be  o ,  Aa ,  oy  ,  and  will  seek  equations  for  the 
increments  in  these  quantities  in  terms  increments  in  the  correspondingly  defined  quantities 
u,Au,uy(=0). 


» 


26 


1 


> 


•  •  • 


•  '• 


I 


In  preparation  for  that,  by  combining  equations  (3-7),  (3-8),  (3-21),  and  (3-22),  the  joint 
response  law  can  be  written  in  the  form 

j  duf+duf  duJN  duff 

au  m  -  =  -  = - =  -  , 

2  sfl  ^2  JIkN 

dbu 1  .  du’  -du’  =  fidu’  -  W'  ♦<<«?)  *  4^-  •  (3-23) 

\J2ks 


(  \2 

.  A  1  A  J 


where  a  superscript  p  denotes  the  contribution  of  "plastic"  or  slipping  joint  behavior.  Similarly, 
by  using  equations  (S-l),  (3-2),  (3-4),  and  (3-22),  the  intact  material  response  can  be  written 


du  I  +  du  J 

du  1  ■  — - - —  -  du  e  +  dulp 

2 

(1  -v)da  -vda 

:  ut  y 


+  wk‘  [aa  -  P(o  +  a  )  -  y] 


(3-24) 


dAu  1  m  duf  -du1  =  dAuIe  +dAulp  =  wi — —  dAo  +  3wkrAa  , 

2  x  £ 

Ie  In  -2 \d5+dOv  . 

du  =  du1  +  dup  =  w - L  +tv?/[aa  -2{kf  -y]  =  0  . 

y  y  y  £  y 

If  the  intact  material  is  in  the  elastic  regime,  equations  (3-24)  still  obtain,  but  with  ?/  -  0  . 

In  the  next  four  subsections  we  will  combine  and  solve  these  equations  for  each  possible 
mode  of  deformation  as  summarized  in  Table  3-3.  Once  the  four  types  of  solutions  are  derived, 
they  are  used  as  the  basis  of  a  computer  program  which  implements  the  logic  embodied  in 
Table  3-4. 


3.5.2.3  Solution  for  mode  1  fEEl  increments.  For  mode  1  (purely  elastic  response), 
equations  (3-23)  with  vanishing  joint  slip  increment  (  dufp  =0  )  and  (3-24)  with  no  intact  rock 
plasticity  (  X1  =  0  )  combine  to  give 


6 


9 


(3-25) 


do  - 


du 

,  dAo  ~ _ 

1  -  V-2V2 

tu 

1 

w. 

E 

v/2*„ 

dov  - 

2 vdo  , 

du!  -  dAa 

y 

S  2 

dAu 


1  +  v 


V2t, 


Note  that  in  mode  1  the  increment  in  joint  shear  displacement  dusJ  is  explicitly  determined  by 
the  joint  constitutive  law  through  the  shear  stiffness  ks. 


» 
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3.5.24  Solution  for  mode  2  fEPI  increments.  For  mode  2  there  is  still  no  intact  rock 
plasticity  (k 1  -  0)  but  now  the  joint  is  slipping  (duJsp  *  0)  .  The  joint  constitutive  law  no 
longer  explicitly  governs  the  incremental  shear  stress-shear  displacement  across  the  joint. 
Instead,  the  Coulomb  friction  law  relates  the  shear  and  normal  stress  across  the  joint: 


fJ  "  (Ts)2  "(HO/J)2  ■  0  ,  (3-26) 

where  the  coefficient  of  friction  is  related  to  the  friction  angle  by  -  tan  <I> .  The  overall 
response  equations  analogous  to  (3-25)  follow  from  (3-23),  (3-24),  and  (3-26): 


► 


> 


do 


du 


1  -  v  -  2v2 
w - 

E 


dAo  -  sgn(Ao)2pcfo 


2 vdo  ,  du 


J 

S 


1  1  +  v  ,  ' 

- \dAu  -w - dAo 


(3-27) 


3.S.2.5  Solution  for  mode  3  rPE'l  increments.  In  mode  3,  the  plastic  proportionality 
constant  k1  becomes  an  additional  unknown.  The  first  three  equations  governing  the  basic 
unknown  increments  are  obtained  as  before  by  combining  (3-23)  and  (3-24),  but  now 
with  k1  *  0  .  An  additional  equation  comes  from  the  condition  of  continued  yielding  (3-6).  In 
matrix  form  the  equations  are 
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0*6 


•  • 


•  • 


P  '• 


1  -v  1 

u>  + 

E 

0 

v 

-w _ 

E 

w^cta  -P(a  +  Oy)  -y] 

do 

du 

0 

w1+v*  1 

0 

3wAa 

d\o 

d&u 

£  v/2*5 

< 

*  -  « 

_  2 
£ 

0 

1 

E 

aoy-2p5-Y 

doy 

0 

2[aa-p(o+oy)-Y] 

3Ao 

2 

ctOy  -2{3o-y 

0 

X7 

0 

(3-28) 

We  have  elected  to  solve  this  set  numerically  during  each  time  step  when  required,  rather  than 
explicitly. 

3 .5.2.6  Solution  for  mode  4  ('PP’l  increments.  In  mode  4  both  the  intact  material  and  the 
joint  are  failing.  The  fundamental  unknowns  are  taken  as  the  three  stress  increments,  X7  ,  and 
the  joint  slip  increment  duJsp .  The  equations  governing  the  stress  increments  follow  as  before 
from  (3-23)  and  (3-24).  The  fourth  equation  is  the  condition  of  continued  intact  material  yielding 
(3-6).  The  fifth  is  the  differential  form  of  the  joint  failure  condition  (3-26),  which  was  already 
used  for  mode  2  and  appears  as  the  second  of  equations  (3-27).  When  assembled  into  matrix 
form  the  equations  are 


1 

0 

v 

~W — 

E 

K-[aa-P(o  +  CT>,)  -y] 

0 

d5 

du 

0 

H.1*’.  1 

0 

3h>Ao 

& 

dho 

dEu 

£  &ks 

_  2 

0 

1 

aay-2Pa-Y 

0 

doy 

■  =  ' 

0 

£ 

£ 

2[aa-p(a  +  ay)-Y] 

3Aa 

___ 

aay  -2^0  “Y 

0 

0 

X7 

0 

-2^sgn(Ac) 

1 

0 

0 

0 

du 

s 

0 

(3-29) 


3.5.3  Analytical  and  Numerical  Results. 

In  addition  to  the  analytic  solution  just  outlined,  complete  solutions  using  explicit  joint 
models  were  submitted  by  CRT,  WA,  and  Itasca.  A  partial  solution  was  completed  by  LLNL. 
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RE/SPEC  encountered  problems  with  the  PRONTO  code — to  be  discussed  later — and  could  not 
generate  a  meaningful  solution.  The  lone  numerical  solution  with  implicit  jointing  was  by  CRT. 
It  gave  results  which  with  one  exception — to  be  discussed  later — were  visually  identical  with 
CRT’s  explicitly  jointed  solution,  so  only  the  latter  will  be  plotted. 

The  zoning  used  by  CRT,  Itasca,  and  LLNL  was  the  simplest  possible  choice — single 
zones  for  each  triangular  block  and  an  explicit  model  of  the  joint.  On  the  other  hand,  WA  broke 
the  structure  into  54  zones  with  six  adjacent  to  the  joint  on  each  side. 

Figures  3-14, 3-15,  and  3-16  contain  the  numerically  predicted  compressibility,  stress  path, 
and  joint  shear  displacement  vs  shear  stress  along  with  the  analytic  solution.  The  figures  show 
that  from  this  perspective  the  numerical  results  that  were  provided  are  fairly  close  to  the  analytic. 
Several  features  of  the  WA  computation  of  net  compressibility  in  Figure  3-14  deserve  mention. 
It  begins  with  the  correct  slope  but  remains  linear  over  most  of  the  duration  of  compressive 
loading.  The  treatment  of  normal  joint  stiffness  was  modified  between  this  problem  and  the 
earlier  ones  such  that  it  was  now  adjustable,  but  still  restricted  to  linear  behavior.  This  solution 
also  contains  much  less  residual  pressure  at  the  end  of  the  calculation,  indicating  that  for  some 
reason  there  was  not  enough  dilatancy  occurring.  Note  also  that  the  LLNL  result  in  Figure  3-16 
is  fairly  accurate,  indicating  that  this  kinematically-based  approach  can  provide  reasonable  results 
for  certain  portions  of  some  problems. 


Volumetric  strain  (percent) 

Figure  3-14.  Overall  compressibility  in  problem  2-S. 
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Figure  3-17  shows  the  deformation  mode  vs  time  for  those  who  provided  this  information. 

It  is  interesting  to  note  that  both  CRT  and  Itasca  predict  a  short  episode  of  mode  4  (PP)  response 

starting  around  /=1.8,  while  the  analytic  solution  does  not.  This  is  where  the  CRT  implicit 

solution  differs  from  explicit.  The  former  stayed  in  mode  4  for  only  one  time  step,  thus  agreeing 

most  closely  with  the  analytic  one.  From  a  practical  standpoint,  in  this  particular  problem,  the 

discrepancy  is  of  little  significance;  the  stresses  and  strains  in  the  four  solutions  are  very  close 

to  one  another.  However,  this  still  raises  some  perplexing  questions.  On  one  hand,  two 

conscientious,  meticulous  calculators  independently  have  reached  the  same  result  with  their 

explicitly  jointed  models;  on  the  other  hand,  the  scheme  detailed  in  Section  3.5.2,  for  numerically  * 

evaluating  the  analytical  solution,  has  been  designed  to  probe  all  possible  evolutionary  paths  from 

each  stress  point  and  to  test  each  one  for  admissibility.  At  every  stage,  this  scheme  has  identified 

a  unique  admissible  stress  increment.  In  particular,  in  the  analytic  solution  (and  effectively  in 

the  CRT  implicit  solution),  the  stress  point  approaches  the  intersection  between  the  two  failure 

surfaces,  moves  right  up  to  it  (by  interpolation),  and  then  moves  directly  across  it,  while  the  two  I 

explicit  numerical  solutions  predict  a  short  traverse  along  the  intersection.  In  the  numerical 

scheme  for  the  analytic  solution,  an  increment  along  the  intersection  (mode  4  or  PP)  was  one  of 

the  four  which  were  explicitly  tested  for  admissibility  (see  the  last  four  rows  of  Table  3-3).  It 
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Figure  3-17.  History  of  deformation  mode  in  problem  2-S. 
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was  found  to  lead  to  a  joint  slip  increment  opposite  the  prevailing  shear  stress  across  the  joint, 
and  so  was  deemed  inadmissible.  It  is  also  interesting  to  note  that  the  elastic  trial  increment 
from  this  point  gave  rise  to  trial  stresses  which  violated  both  failure  conditions  (and  was  ruled 
out  for  that  reason).  This  is  not  the  same  as  the  elastic  trial  used  by  the  numerical  solutions4. 
Nevertheless,  one  might  speculate  that  the  numerically  derived  trial  stresses  in  the  explicit  models 
also  violated  both  failure  conditions  and  that  this  somehow  initiated  the  excursion  along  the 
intersection  of  the  failure  surfaces.  It  may  then  have  taken  these  methods  a  number  of  time  steps 
to  recognize  and  correct  for  joint  slippage  opposite  the  shear  stress  across  the  joint. 

3.5.4  Why  PRONTO  Was  Not  Applicable  to  Problem  2-S. 

RE/SPEC  elected  always  to  treat  the  elastic  properties  of  joints  implicitly,  i.e.,  by 
smearing  them  out  and  distributing  them  in  the  surrounding  intact  material.  Only  the  Coulomb 
slip  is  represented  explicitly.  This  may  be  reasonable  if  there  are  many  such  joints,  but  with  just 
a  single  one  it  leads  to  difficulties  when  the  joint  is  inclined.  The  problem  stems  from  the  fact 
that  the  RE/SPEC  approach  makes  the  surroundings  behave  like  a  transversely  isotropic  material, 
with  the  primary  material  axis  aligned  normal  to  the  joints).  In  problem  2-S  this  axis  is  45 
degrees  from  the  geometric  axes.  To  see  the  effect  of  this  misalignment,  consider  the  situation 
before  any  slippage  occurs.  The  two  triangular  blocks  are  acting  as  a  single,  continuous  square 
one.  Because  of  the  anisotropy,  the  specified  boundary  conditions  (constant  normal 
displacements,  vanishing  shear  traction)  are  not  consistent  with  a  homogeneous  stress  state  within 
the  block.  In  other  words,  with  no  shear  traction  the  block  would  suffer  a  shear  strain  and  would 
require  linearly  varying  normal  boundary  displacements  to  maintain  homogeneity;  conversely, 
a  particular  value  of  shear  traction  could  be  applied  which  would  cause  shear  strain  to  vanish. 
Therefore,  when  RE/SPEC  set  up  and  ran  this  problem  with  the  specified  boundary  conditions, 
the  stress  state  was  nonuniform  from  the  very  beginning.  Later  the  situation  got  even  worse, 
because  the  joint  failure  condition  was  not  met  simultaneously  all  along  the  joint,  not  even 
approximately.  Because  the  nonuniformities  were  traceable  to  particular  features  of  the  RE/SPEC 
approach  in  a  way  that  was  not  anticipated  when  the  problem  was  defined,  it  was  decided  not 
to  pursue  the  solution  to  its  conclusion  with  this  method. 


4For  this  problem,  the  trial  elastic  stresses  used  in  the  numerical  solutions  come  from  strains 
within  each  element  which  depend  in  part  on  internal  (non-forced)  nodal  displacements,  projected 
from  the  previous  time  step  by  forward  differences  based  on  the  momentum  balace  equations. 
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SECTION  4 


PROBLEM  STATEMENTS  AND  RESULTS  FOR  TWO-DIMENSIONAL  PROBLEMS 


4.1  PRELIMINARY  DISCUSSION. 

The  problems  treated  so  far  could  be  classified  as  both  static  and  "zero-dimensional," 
since  their  underlying  solutions  have  neither  essential  temporal  nor  large-scale  spatial  variation. 
As  we  have  seen,  this  simplicity  makes  analytic  solutions  possible,  and  those  have  been  used  as 
ground  truth  for  assessing  the  numerical  results.  In  contrast,  problems  3  and  4  have  fields  which 
vary  both  spatially  and  temporally.  Analytic  solutions  are  not  generally  feasible,  and  the  proper 
procedure  for  assessing  the  accuracy  of  numerical  solutions  is  much  less  obvious.  In  this  study, 
we  have  tried  to  use  all  of  the  tools  available  to  evaluate  the  numerical  solutions,  but  none  of 
them  is  as  definitive  as  a  full  analytic  solution  would  be.  So  we  check  for  internal  consistency, 
compatibility  with  the  basic  understanding  of  wave  propagation  processes,  robustness  of  the 
method  as  revealed  in  the  zero-dimensional  problems,  and  comparison  with  the  other  solutions. 
We  can  perform  limited  analytical  solutions,  for  example  taking  numerically  computed 
macroscopic  strains  as  input  and  deriving  intact  rock  response  analytically.  In  the  end  we  will 
form  some  judgments  about  the  relative  credibility  of  the  various  solutions  to  these  particular 
problems.  But  such  judgments  cannot  ever  be  definitive  for  problems  with  no  complete  analytic 
solution.5 

4.2  PROBLEM  3:  DIVERGENT  WAVE  PROPAGATION  THROUGH  A  JOINTED  ROCK 

ISLAND. 

4.2.1  Statement  of  the  Problem. 

This  problem  concerns  deformations  of  a  wedge-shaped  section  of  an  annulus  in  plane 
strain,  as  shown  in  Figure  4-1.  The  entire  region  contains  vertically  and  horizontally  jointed  rock 
as  specified  in  Tables  2-1  and  2-2.  The  top  edge  (the  inner  arc)  is  loaded  with  the  pressure  pulse 
shown  in  the  figure,  while  shear  tractions  are  zero.  (The  pulse  approximates  that  at  a  500-m 
range  from  1/3  MT  of  coupled  nuclear  energy).  The  left  and  right  sides  have  roller  boundaries, 
making  the  left  side  a  plane  of  symmetry  (the  right  side  is  not,  because  the  effective  anisotropy 
due  to  jointing  makes  the  material  unsymmetric  about  that  plane).  The  lower  edge  (the  outer  arc) 
is  a  transmitting  boundary.  Most  of  the  region  is  to  be  modelled  implicitly,  except  for  a 
rectangular  region  extending  2.5  tunnel  diameters  (12.5  m)  in  all  directions  from  the  on-axis  point 
at  R=500  m.  This  region  is  to  be  modelled  explicitly,  in  anticipation  of  the  tunnel  situated  here 
in  problem  4.  The  problem  is  to  find  stresses,  velocities,  and  strains  throughout  the  region. 


5Some  individuals  believe  that  numerical  methods  can  be  "validated"  by  comparison  with 
experiments.  Whether  or  not  this  is  true  is  mainly  a  matter  of  definition.  But  in  this  study,  using 
data  is  not  an  option.  All  the  material  models  are  highly  idealized  and  not  intended  to  accurately 
represent  real  material  behavior. 


.Source 


Implicitly  Jointed 
Rock 


Transmitting  Boundary 


Time  (ms) 


Figure  4-1.  Geometry  and  loading  in  problem  3. 


When  this  problem  was  originally  set  up,  a  gravitational  prestress  was  proposed.  This 
seemingly  innocuous  requirement  led  to  both  philosophical  disagreements  and  computational 
difficulties,  and  was  eventually  discarded.  In  principal,  given  a  material  model  and  some 
assumption  about  the  lateral  lithostatic  stress,  it  is  straightforward  to  define  the  gravitationally 
induced  stresses  at  any  depth  in  a  half-space,  then  translate  those  to  a  set  of  boundary  tractions, 
body  forces,  and  stresses  in  the  wedge-shaped  region  of  problem  3.  Modelers  should  have  been 
able  to  impose  these  loads  as  initial  conditions,  then  proceed  forward  with  the  active  loads  and 
boundary  conditions  applied  relative  to  the  initial  state. 

In  practice,  most  of  the  modelers  were  unable  to  proceed  as  above  without  either 
extensive  code  modifications  or  separate  off-line  calculations.  Some  could  not  impose  body 
forces,  some  could  not  impose  the  necessary  shear  tractions  on  those  boundary  segments  that 
were  neither  vertical  nor  horizontal,  and  some  could  not  switch  the  character  of  the  boundary 
conditions  between  the  initialization  phase  and  the  dynamic  loading.  A  compromise  approach 
was  to  ignore  body  forces  and  boundary  shears,  and  impose  only  a  normal  traction  on  the  upper 
boundary  arc  as  the  initial  condition.  But  even  if  the  lower  boundary  could  have  been  made  to 
switch  character  (statically  a  roller,  dynamically  transmitting),  the  vanishing  shears  on  the  slanted 
right-hand  edge  and  upper  arc  would  have  led  to  an  initial  stress  state  unlike  the  intended  one. 


Furthermore,  even  if  a  reasonable  initial  stress  state  could  have  been  defined,  the  introduction  of 
a  tunnel  in  problem  4  would  have  raised  a  whole  new  set  of  issues,  differences  of  opinion,  and 
difficulties.  For  example,  one  point  of  contention  would  have  been  whether  to  impose  the 
prestress  after  the  tunnel  is  emplaced,  or  to  computationally  "excavate"  the  tunnel  in  an  already 
prestressed  rock  mass.  In  light  of  all  the  potential  pitfalls,  the  limited  insight  into  real 
computational  differences  that  would  have  been  gained  compared  with  the  effort  that  would  have 
been  required,  and  the  small  anticipated  impact  on  the  results  even  if  the  loads  were  imposed 
correctly,  the  gravitational  prestress  was  abandoned. 

4.2.2  Physical  Effects  Observed  in  the  Numerical  Solutions. 

The  gross  features  of  the  numerical  solutions  are  quite  consistent  with  each  other  and  with 
expectations.  Above  all,  there  is  very  little  variation  in  the  stress  and  deformation  pulses  over 
the  computational  grid,  except  for  temporal  offsets  which  are  fully  consistent  with  anticipated 
wave  propagation  speeds.  The  on-axis  vertical  stress  and  velocity  at  the  center  of  the  grid  are 
shown  in  Figures  4-2  and  4-3.  Note  the  small  attenuation  in  peak  stress  that  has  occurred  from 
the  top  edge,  where  the  applied  peak  was  100  MPa,  to  the  center,  where  peaks  of  99.1,  94.1, 
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Figure  4-2.  Radial  stress  at  tunnel  location  in  problem  3. 
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Figure  4-3.  Radial  velocity  at  tunnel  location  in  problem  3. 

94.7,  and  100.8  MPa  were  reported  by  CUT,  WA,  Itasca,  and  RE/SPEC.  Attenuation  can  be 
traced  either  to  geometric  divergence  or  material  dissipation.  The  pure  geometric  divergence  of 
the  region  over  this  range  is  slight:  R  _1/2  is  only  5  percent  smaller  at  the  center  than  the  top. 
Furthermore,  the  effective  anisotropy  induced  by  the  joints  might  be  focussing  energy  in  the 
downward  direction. 

Note  that  the  LLNL  approach  produced  a  fairly  accurate  velocity  pulse.  This  is  probably 
a  reflection  of  the  fidelity  of  the  overali  uniaxial  compressibility  as  lumped  into  the  joints. 

Another  indicator  of  divergence  is  the  post-peak  shape  of  the  strain  path.  Figure  4-4 
shows  calculated  strain  paths  at  R=500  m,  0=2 ,  in  implicitly  jointed  rock.  They  present  a 
macroscopic  view  of  the  deformation  of  a  homogeneous  effective  medium  with  the  properties  of 
intact  rock  and  joints  combined  and  smeared  out.  In  this  medium,  as  in  any  homogeneous 
medium  under  rapidly  rising  but  slowly  decaying  divergent  dynamic  loading,  the  strain  path  is 
almost  purely  uniaxial  up  to  the  initial  peak  (i.e.,  very  little  hoop  strain),  but  on  unloading  the 
hoop  strain  becomes  increasingly  negative  (tensile)  as  outward  radial  displacement  accumulates. 
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Figure  4-4.  Strain  path  at  500-m  range,  2  degrees  off  axis,  in  problem  3. 


Because  the  medium  is  a  composite  (of  joints  and  intact  rock),  the  foregoing  observations 
do  not  fully  reveal  the  material’s  deformation  under  divergent  wave  loading.  In  particular,  while 
the  macroscopic  response  is  mainly  uniaxial  strain  on  loading,  the  intact  rock  constituent 
undergoes  substantial  hoop  expansion  at  the  same  time.  This  is  illustrated  in  Figure  4-5,  which 
shows  Itasca’s  calculated  on-axis  strain  paths  both  above  and  below  the  upstream  and 
downstream  edges  of  the  explicitly  modelled  region.  In  the  implicitly  jointed  material  the  first 
leg  is  uniaxial  strain  compression;  the  second,  simultaneous  hoop  and  radial  expansion,  and  the 
third  (with  some  imagination)  pure  radial  expansion.  On  the  other  hand,  in  the  intact  rock  most 
of  the  hoop  expansion  occurs  during  the  first  leg  along  with  the  radial  compression.  This  effect 
is  exaggerated  by  the  particular  constraint  conditions  of  this  planar  model,  and  may  not  be 
generally  representative  of  three-dimensional  divergent  wave  propagation.  Here,  through  a 
Poisson  effect,  the  out-of-plane,  plane  strain  constraint  makes  the  intact  rock  behave  more  stiffly 
in  the  plane.  When  compressed  radially  (vertically)  it  develops  compressive  hoop  (horizontal) 
stresses  which  are  larger  than  they  would  be  in  the  absence  of  the  out-of-plane  constraint.  The 
larger  compressive  stresses  are  transmitted  directly  into  the  vertical  joints.  Because  the 
constrained  intact  rock  is  stiffer  than  the  joints,  the  joints  compress  and  the  intact  rock  expands 
horizontally. 
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Figure  4-5.  On-axis  strain  paths  in  implicitly  jointed  and  intact  rock,  analytic  solution  and 
Itasca  calculation,  problem  3. 


Figure  4-5  also  shows  an  analytic  representation  of  this  effect.  The  light  solid  curve 
labelled  "input  for  analytic  solution"  is  simply  an  idealization  of  the  two  Itasca  strain  paths  in 
implicitly  jointed  rock.  This  path  was  used  as  input  into  a  calculation  using  the  machinery 
developed  for  problem  1-IM,  i.e.,  a  homogeneous  block  of  implicitly  jointed  material  was 
homogeneously  strained  according  to  the  prescribed  path.  The  heavy  solid  line  is  the 
corresponding  strain  path  in  intact  rock.  Note  that  it  agrees  well  with  the  numerical  results  in 
both  magnitude  and  shape.  The  kink  in  the  loading  leg  signals  the  onset  of  plastic  flow.  The 
fact  that  the  two  kinds  of  paths  end  up  at  about  the  same  point  is  related  to  the  fact  that  during 
this  deformation  the  joint  behaves  elastically,  so  when  the  stresses  become  very  small  the  joints 
return  to  near  their  original  dimensions  and  the  residual  strains  are  mostly  confined  to  the  intact 
rock. 


Returning  to  the  matter  of  dissipation,  in  this  problem  it  could  arise  either  from  intact 
rock  plasticity  or  joint  slip.  The  kinematic  constraints  and  loading  virtually  preclude  joint  slip. 
Plasticity  begins  part  way  through  the  loading  phase,  which  macroscopicaily  is  essentially 
uniaxial  strain.  On  unloading,  the  small  geometric  divergence  does  increase  the  deviatoric  stress 


over  that  which  would  occur  in  pure  uniaxial  strain,  and  this  enhances  the  plasticity.  To  illustrate 
this,  Figure  4-6  shows  the  stress  path  at  a  typical  on-axis  point  from  all  the  numerical 
calculations,  compared  with  the  analytic  result  for  uniaxial  strain  up  to  about  the  same  peak 
strain.  (The  latter  comes  from  the  general  solution  for  problem  2.)  For  both  strain  paths 
plasticity  begins  at  a  pressure  of  about  15  MPa  and  continues  up  to  peak  stress.  But  on 
unloading,  the  deviator  stress  is  much  greater  at  the  same  pressure  in  the  divergent 

problem  than  the  uniaxial  one.  The  divergence,  though  small,  is  sufficient  to  cause  the  unloading 
to  be  predominantly  plastic,  following  along  the  failure  surface,  while  by  comparison  it  remains 
elastic  in  the  uniaxial  case. 
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Figure  4-6.  Stress  path  at  487-m  range  on  axis  in  implicitly  jointed  rock  in  problem  3. 


With  all  this  plasticity,  why  then  isn’t  there  more  attenuation  due  to  dissipation?  The 
main  reason  is  dilatancy,  which  tends  to  offset  the  deviatoric  dissipative  losses  with  volumetric 
expansion  against  a  positive  pressure.  This  is  illustrated  in  Figure  4-7,  which  gives  the  pressure- 
volume  paths  computed  in  adjoining  zones  in  the  implicitly  and  explicitly  jointed  rock.  Dilatancy 
causes  a  larger  final  volume  than  initial,  corresponding  to  negative  volumetric  dissipation  (the 
area  between  the  loading  and  unloading  curves).  Net  volumetric  work  is  being  done  by  the 
material.  The  difference  between  the  two  sets  of  curves  in  Figure  4-7  is  that  the  one  in  the 
explicitly  jointed  region  represents  local  volumetric  strain  in  the  intact  rock  only,  white  the  other 
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Figure  4-7.  Compressibility  at  487-m  range  on  axis  in  problem  3. 


set  shows  overall  volumetric  strain  including  the  portion  due  to  the  joints.  (There  is  only  one 
value  of  pressure;  it  is  the  same  in  intact  rock  as  in  implicitly  jointed  rock.)  The  volumetric 
strain  due  to  the  joints  is  entirely  elastic,  positive,  and  recoverable  (compressive),  while  the  intact 
material  contributes  a  positive  elastic  portion  and — after  the  onset  of  plasticity  at  about  15  MPa 
pressure — a  negative  (expansive)  plastic  dilatant  portion.  Under  the  constraint  conditions  here, 
at  low  pressure  the  elastic  volumetric  strain  in  intact  rock  is  about  one  third  of  the  total  elastic 
strain.  This  fraction  increases  as  the  joints  stiffen  with  pressure  (in  the  limit  of  very  large 
pressure  the  joints  would  become  rigid  and  all  the  elastic  strain  would  be  in  the  intact  rock).  The 
dilatancy  accumulates  continuously  during  both  the  post-yield  portion  of  the  loading  phase  and 
the  plastic  portions  of  the  unloading.  Referring  to  Figure  4-7,  during  loading  the  dilatant  volume 
strain  is  the  horizontal  distance  between  the  actual  curve  and  the  extrapolation  of  the  initial 
(elastic)  part  to  higher  pressures.  During  unloading,  the  additional  dilatant  volume  strain  is  the 
horizontal  distance  between  the  unloading  and  loading  curves.  At  the  same  piessure,  these 
horizontal  distances  should  be  the  same  in  the  implicit  curve  as  the  intact  one,  since  there  is  no 
plasticity  in  the  joints.  Consequently,  the  areas  representing  work  done  by  the  material  will  be 
the  same  regardless  of  which  type  of  pressure-volume  curve  is  used. 
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4.2.3  Numerical  Issues. 


This  problem  was  specified  to  have  a  nonreflecting  boundary  at  the  lower  edge  of  the 
grid.  It  is  well  known  that  no  all-purpose  absorbing  boundary  is  possible.  However,  in  this  * 

problem  practically  all  of  the  motion  is  normal  to  the  boundary  and  the  material  is  almost  elastic, 
so  in  principle  it  is  possible  to  meet  the  specified  condition  very  accurately.  One  obvious  way 
to  evaluate  the  absorbing  boundaries  is  to  examine  stress  pulses  at  various  locations  for  reflected 
signals.  Figure  4-8  shows  radiai  stress  pulses  at  three  ranges  for  each  calculator  who  submitted 
the  data.  The  data  from  RE/SPEC  show  clear  evidence  of  a  fairly  strong  reflection,  and  a  little  > 

analysis  of  the  timings  shows  that  it  has  to  come  from  the  bottom  edge  of  the  grid.  The 
explanation  is  that  the  PRONTO  formulation  requires  compressional  and  shear  wave  speeds,  and 
computes  them  without  accounting  for  the  anisotropy  introduced  into  the  implicitly  modelled 
surroundings  by  RE/SPEC’s  "compliant  joint"  model.  Thus  with  incorrect  wave  speeds  a 
spurious  reflection  is  generated.  RE/SPEC  did  not  directly  remedy  this  error  during  the  course 
of  this  study,  although  this  problem  was  rerun  with  an  extended  grid  such  that  no  bottom 
boundary  reflection  reached  the  region  of  interest  during  the  time  of  interest. 


Figure  4-8.  Radial  stress  pulse  shapes  at  480,  500,  and  520-m  ranges  in  problem  3. 
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A  related  issue  concerns  the  continuity  between  the  implicitly  and  explicitly  modelled 
regions.  In  principle  the  boundaries  should  be  invisible  to  waves  which  are  long  compared  to 
joint  spacing.  In  this  problem  the  rise  time  is  10  ms  and  the  loading  wave  speed  is  about 
2.5  m/ms,  so  the  stress  rise — the  shortest  feature  of  the  wave — is  spread  out  over  about  25  m. 
Since  the  joints  are  only  1  m  apart  one  would  not  expect  to  see  much  effect  of  the 
implicit/expiicit  boundary.  In  Figure  4-8  there  is  no  obvious  sign  of  reflections  from  these 
boundaries.  Another  diagnostic  of  possible  material  discontinuity  is  the  difference  in  computed 
horizontal  stress  pulses  on  either  side  of  the  boundary;  they  should  be  equal.  These  traces  at  the 
downstream,  on-axis  discontinuity  are  shown  in  Figure  4-9.  None  of  the  calculations  gives  a 
perfect  match,  but  CRT  and  RE/SPEC  are  quite  close,  while  WA  and  Itasca  differ  by  20  percent 
or  more  over  most  of  the  pulse.  The  discrepancy  may  be  related  to  the  fact  that  both  WA’s  and 
Itasca’s  implicit  models  are  isotropic  (Appendices  D  and  E),  while  as  noted  above  the  overall 
behavior  of  the  jointed  rock  mass  is  orthotropic.  Note,  however,  that  by  adjusting  Poisson’s  ratio 
even  an  isotropic  material  could  be  made  to  give  nearly  the  correct  coupling  between  horizontal 
and  vertical  stresses  and  strains,  and  therefore  to  provide  the  same  horizontal  stresses  in  the 
implicitly  and  explicitly  jointed  regions  when  the  wave  traverses  that  boundary.  Apparently  WA 
and  Itasca  used  some  other  criterion  for  selecting  the  material  properties  in  the  implicitly  jointed 
region.  The  orthotropic  elastic  idealization  is  discussed  further  in  Section  4.3.3. 


Figure  4-9.  Hoop  stress  at  512-m  range  on  axis,  implicit  and  explicit  region,  problem  3. 


4.3  PROBLEM  4:  LINED  TUNNEL  IN  A  JOINTED  ROCK  ISLAND  SUBJECTED  TO 

DIVERGENT  WAVE  PROPAGATION 

4.3.1  Statement  of  the  Problem. 

The  geometry  and  loading  in  this  problem  are  shown  in  Figure  4-10.  They  are  precisely 
the  same  as  problem  3,  but  now  there  is  a  lined  tunnel  in  the  center  of  the  rock  island.  The 
properties  of  the  liner  have  been  listed  in  Tables  2-1  and  2-2.  Calculators  were  asked  to  provide 
many  diagnostics  of  tunnel  and  surrounding  medium  response,  some  of  which  will  be  discussed 
here. 

The  LLNL  methodology  has  no  intrinsic  mechanism  for  representing  Poisson  effects,  i.e., 
horizontal  stresses  caused  by  vertical  strains.  Therefore,  in  some  cases,  two  different  solutions 
are  presented,  one  with  no  horizontal  stress,  and  a  second  with  an  externally  imposed,  constant, 
horizontal  stress  of  33  MPa.  Unless  otherwise  noted,  the  LLNL  results  are  for  the  nonvanishing 
horizontal  stress. 


Figure  4-10.  Geometry  and  loading  in  problem  4. 


4.3.2  Physical  Effects  Observed  in  the  Numerical  Solutions. 

To  begin  this  discussion,  some  general  observations  will  be  made.  The  incident  pulse  was 
analyzed  in  problem  3  and  is  shown  in  Figures  4-2,  4-3,  and  4-8.  It  propagates  at  about 
2.5  m/ms  with  little  change  in  shape.  Thus  the  leading  edge  of  the  pulse  arrives  at  the  tunnel 
location  about  20  ms  after  application  to  the  rock  island  boundary,  and  the  peak,  about  10  ms 
later. 

It  has  been  noted  both  in  experiments  and  calculations  that  tunnel  deformations  are 
generally  more  severe  under  divergent  wave  loading  than  uniaxial  loading  to  the  same  peak 
stress.  Conventional  wisdom  holds  that  in  divergent  flow  the  post-peak  radially  outward  motions 
cause  rapid  circumferential  unloading,  "loss  of  confinement,"  and  consequent  reduced  strength 
around  the  tunnel,  compared  with  unidirectional  flow.  We  have  already  seen  in  problem  3  that 
the  macroscopic  free  field  strain  path  in  this  geometry  does  include  post  peak  hoop  expansion 
(Figures  4-4,  4-5),  and  that  the  stress  path  remains  either  on  the  failure  surface  on  unloading  or 
much  closer  to  it  than  in  uniaxial  strain  (Figure  4-6).  However,  according  to  Figures  4-11  and 
4-12,  which  show  histories  of  calculated  crown/invert  and  springline  closures,  all  participants 
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Figure  4-11.  Crown-invert  closure  in  problem  4. 
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Figure  4-12.  Springline  closure  in  problem  4. 
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except  RE/SPEC  predict  that  before  35  ms  have  elapsed,  the  tunnel  deformation  has  essentially 
reached  its  transient  maximum.  At  this  time,  only  about  20  percent  of  the  ultimate  outward 
displacement  and  hoop  strain  has  accumulated,  and  the  free  field  stress  at  the  tunnel  range  has 
only  dropped  by  about  5  percent  from  its  peak.  Moreover,  results  from  problem  3  indicate  that 
after  the  peak,  the  stress  point  is  often  slightly  below  the  failure  surface  (Figure  4-6),  not 
precisely  on  it  as  the  "loss-of-confinement"  scenario  would  seem  to  imply.  Thus,  while  some 
aspects  of  the  calculated  results  appear  to  conform  with  conventional  wisdom,  other  details 
suggest  that  the  situation  is  more  complex. 

One  crucial  feature  which  bears  on  the  relationship  between  the  free  field  and  the  tunnel 
deformation  is  the  stress  and  strain  concentration  caused  by  the  presence  of  the  tunnel  itself.  For 
example,  if  the  surroundings  were  infinite,  isotropic,  and  elastic  and  the  free  field  stress  state 
static  and  unidirectional,  then  it  is  well  known  that  the  tangential  stress  at  the  springline  of  a 
circular  tunnel  is  three  times  the  free  field  stress.  While  plasticity  and  joint  failure  place  limits 
on  the  achievable  stress  concentration  in  this  problem,  it  is  clear  that  the  free  field  stress  state 
need  not  be  precisely  on  the  failure  surface  in  order  for  the  material  near  the  tunnel  to  in  a  plastic 
state. 
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Figures  4-13  to  4-16  contain  a  more  detailed  picture  of  the  predicted  stresses  around  the 
tunnel.  They  show  the  radial  and  tangential  (with  respect  to  the  tunnel)  stress  distributions  near 
the  tunnel  at  two  different  times,  free  field  peak  stress  arrival  (about  30  ms)  and  end  of  positive 
phase  (about  120  ms).  Four  of  the  calculators  (CRT,  RE/SPEC,  and  Itasca)  provided  results  in 
this  form.  The  curves  labelled  "elastic”  in  the  first  two  figures  will  be  discussed  later. 

First  consider  the  radial  stresses  (with  respect  to  tunnel  centerline)  at  peak  stress  arrival 
(Figure  4-13).  According  to  the  results  of  problem  3  the  free  field  vertical  stress  at  the  tunnel 
range  is  about  100  MPa  and  the  horizontal  stress  about  30  MPa.  At  this  time  the  free  field  I 

stresses  will  decrease  with  either  increasing  or  decreasing  range  from  the  source  point,  by  about 
40  or  5  percent  respectively  over  the  12.5  m  (2.5  diameters)  plotted.  While  the  calculated  radial 
stresses  12.5  m  above  the  tunnel  do  exceed  those  an  equal  distance  below  (as  expected),  neither 
level  is  as  high  as  the  corresponding  free  field  value.  This  indicates  that  the  relief  provided  by 
the  tunnel  wall  has  an  effect  at  least  this  far  out  in  the  vertical  direction.  Along  a  horizontal  ^ 

radial,  however,  the  calculated  radial  stresses  approach  the  free  field  horizontal  stress  of  30  MPa 
much  more  quickly. 

Note  in  Figure  4-13  that  radial  stresses  along  all  radials  approach  about  8  MPa  at  the 
tunnel  wall.  This  is  the  pressure  exerted  by  the  yielding  perfectly  plastic  liner  in  the  "breathing" 
mode.  * 

Figure  4-14  clearly  shows  the  concentration  of  tangential  stress  along  the  horizontal  radial. 

Peak  reported  values  range  between  120  and  160  MPa,  compared  to  the  free  field  stress  in  the 

same  direction  of  about  100  MPa.  However,  these  peaks  are  not  as  high  as  an  isotropic  elastic 

analysis  would  predict.  The  reason,  as  noted  above,  is  that  the  plasticity  limits  the  level  of  » 

achievable  stress  concentration,  particularly  in  the  presence  of  the  relief  afforded  by  the  tunnel 

wall. 

The  "sawtooth"  variation  of  horizontal  stress  along  the  vertical  radials  in  Figure  4-14  is 
due  to  target  point  locations  alternating  between  intact  rock  and  joints.  There  is  no  requirement  ^ 

that  horizontal  stresses  in  horizontal  joints  be  continuous  with  those  in  adjoining  intact  rock.  It 
is  comforting  to  note  that  most  calculations  approximately  agree  on  the  shapes  and  magnitudes 
of  these  curves. 

Thus  many  features  of  the  stress  states  at  30  ms  can  easily  be  explained  based  on 
intuitive,  static  analysis,  and  more  will  be  explained  later  based  on  an  orthotropic  elastic  analysis.  * 

In  contrast,  the  picture  at  the  end  of  the  positive  phase  (Figures  4-15  and  4-16)  is  not  nearly  as 
clear.  The  only  stresses  present  are  either  dynamic  or  residual,  and  the  levels  are  much  lower 
on  the  average.  One  reason  for  RE/SPEC’s  larger  values  compared  with  the  others  is  that  they 
plotted  results  at  the  end  of  their  free-field  positive  velocity  phase,  which  came  earlier  than  the 
others  (see  Figure  4-3),  while  the  free-field  stress  was  still  moderately  compressive  (see  > 

Figure  4-2). 

Figures  4-17  and  4-18  show  the  slip  along  the  first  four  horizontal  joints  above  the  tunnel 
centerline  at  two  times,  peak  free-field  stress  arrival  and  end  of  positive  phase.  (All  calculators 
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Figure  4-13.  Radial  stress  field  near  the  tunnel  at  peak  free  field  stress  arrival  time  (about 
30  ms). 
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tial  stress  field  near  the  tunnel  at  peak  free  field  stress  arrival  time  (about 
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Figure  4-16.  Tangential  stress  field  near  the  tunnel  at  the  end  of  the  free-field  positive  phase 
(about  120  ms). 
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Figure  4-18.  Slip  along  horizontal  joints  at  the  end  of  the  free-field  positive  phase. 
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predicted  nearly  equal  and  opposite  slip  on  the  corresponding  joint  below  the  centerline  at  the 
same  time.)  Virtually  all  the  slips  are  in  the  sense  that  would  shorten  the  springline.  The 
greatest  amount  of  slip  occurs  along  the  joint  tangent  to  the  top  of  the  tunnel  (joint  2)  and  the 
one  below  that  (joint  3).  The  slip  at  the  later  time  exceeds  that  at  the  earlier  one  by  an  average 
factor  of  about  four.  At  the  earlier  time,  CRT  appears  to  have  some  non-slipping  nodes 
surrounded  by  slipping  ones,  which  runs  contrary  to  intuition.  Itasca  and  WA  predict  slips  that 
average  2  to  3  times  the  others,  but  CRT  predicts  a  large  slip  on  joint  3  at  the  end  of  the  positive 
phase  at  the  tunnel  wall.  Overall  there  is  very  little  consistency  in  these  results  among  the 
calculators,  making  them  very  difficult  to  interpret. 

One  question  that  could  uikcd  is  "How  much  springline  closure  can  be  traced  directly 
to  joint  slip?"  Table  4-1  addresses  this  question.  The  table  simply  sums  the  slips  estimated  at 
the  tunnel  wall  at  the  end  of  the  positive  phase  for  the  two  joints  nearest  the  centerline  (joints 
3  and  4).  These  two  joints  could  be  said  to  contribute  directly  to  springline  closure  since  they 
are  the  ones  that  intersect  the  upper  half  of  the  tunnel  boundary.  Again  there  is  little  consistency. 
CRT  and  Itasca’s  joints  contribute  to  order  unity,  while  RE/SPEC’s  barely  contribute  at  all. 
However,  this  ranking  is  not  representative  of  the  average  slips,  where  CRT  and  RE/SPEC  are 
about  equal,  but — contrary  to  intuition — RE/SPEC’s  predictions  appear  to  go  to  zero  at  the  tunnel 
wall. 


Table  4-1.  Joint  slip  contribution  to  springline  closure. 


CRT 

RE/SPEC 

Itasca 

WA 

Joint  3  slip  at  tunnel  edge  (m) 

0.0904 

0.0060 

0.0278 

0.0169 

Joint  4  slip  at  tunnel  edge  (m) 

0.0023 

0.0004 

0.0152 

0.0047 

Sum  of  joint  3  and  4  slips  (m) 

0.0927 

0.0064 

0.0430 

0.0216 

Half  of  springline  closure  (m) 

0.138 

0.325 

0.165 

0.145 

Fraction  of  closure  due  to  slip 

0.67 

0.02 

0.26 

0.15 

Numerically  predicted  final  deformed  tunnel  shapes  are  shown  in  Figure  4-19 
(Displacements  have  been  exaggerated  by  a  factor  of  10).  According  to  this  figure  and 
Figures  4-11  and  12,  all  calculations  except  LLNL’s  show  reverse  ovalling,  i.e.,  the  springline 
closure  is  positive.  Most  calculations  predict  severe,  localized  deformation  of  the  block  of  rock 
bordering  the  springline.  This  appears  to  be  due  to  a  combination  of  three  factors:  (1)  Highly 
concentrated  vertical  stresses  at  the  springline  combine  with  (2)  low  lateral  confinement  due  to 
the  presence  of  the  tunnel  wall  (and  not  necessarily  to  wavefront  divergence)  to  create  a  failed, 
vertically  compressed  stress  state.  (3)  Dilatancy  then  inhibits  net  compaction,  exaggerating  the 
amount  of  "extrusion"  of  intact  rock  into  the  tunnel.  Note  that  since  the  LLNL  model  did  not 
represent  the  rock  as  a  deformable  continuum,  it  could  not  have  accurately  represented  the  effect 
of  dilatancy,  nor  could  it  have  correctly  predicted  the  stress  concentration  near  the  tunnel  wall. 


r 


4.3.3  Orthotropic  Elastic  Predictions  of  Stresses. 

At  a  propagation  speed  of  about  2.5  m/ms,  the  rising  portion  of  the  incident  pulse  extends 
over  about  25  m  and  the  decaying  portion  over  more  than  200  m.  Both  of  these  distances  are 
fairly  large  compared  with  the  tunnel  diameter  of  5  m,  so  there  is  a  possibility  that  some  aspects 
of  the  response  will  be  quasi-static.  Moreover,  as  noted  before,  the  dilatancy  causes  the  plastic 
response  to  be  "less  plastic,"  and  in  particular  there  is  less  softening  on  initial  loading  above  the 
yield  point  than  there  would  otherwise  be.  So  we  might  expect  that  during  initial  loading,  a  two- 
dimensional,  static,  orthotropic,  elastic  solution  could  shed  some  light  on  the  stress  distributions 
around  the  tunnel.  To  this  end,  in  this  section  we  will  construct  such  a  solution  under  the 
loading  conditions  in  force  when  the  peak  incident  stress  just  reaches  the  tunnel  center. 

The  derivation  of  the  orthotropic  elastic  solution  follows  Lekhnitski  (1966).  Ke  explains 
the  general  approach  to  two-dimensional  problems  using  complex  potentials,  and  provides  some 
examples  for  infinite  regions  with  holes,  but  does  not  specifically  work  out  the  full  stress  fields 
for  the  fundamental  problem  of  interest  here,  viz.,  an  infinite  region  with  a  hole  and  with  uniaxial 
stress  at  infinity.  So  we  will  fill  in  the  details  needed  to  derive  the  fundamental  solution,  and 
show  how  this  solution  can  be  rotated  and  superposed  to  build  up  one  containing  the  effects  of 
vertical  and  horizontal  free-field  stress  and  internal  pressure  to  represent  the  effect  of  the  liner. 

This  first  matter  to  be  addressed  is  the  choice  of  orthotropic  elastic  constants.  In  this 
section  only,  let  x,  y,  and  z  represent  the  horizontal,  vertical,  and  out-of-plane  directions 
respectively.  The  elastic  relation  among  macroscopic  normal  stresses  and  strains  in  implicitly 
jointed  material  is  similar  to  (3-18): 
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We  are  interested  in  plane  strain  in  the  z-direction,  so  we  set  ez  =  0  .  The  resulting  expressions 
for  normal  strains,  when  compared  with  the  standard  forms 
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For  shear  components,  we  write  the  engineering  shear  strain  y  as  a  sum  of  two  parts,  the 
first  x^G  due  to  intact  rock  shearing,  and  the  second  due  to  joint  elastic  shear.  For  the  latter, 
note  that  a  square  block  of  intact  rock  of  any  size  w,  bounded  by  joints,  will  suffer  relative 
tangential  displacements  xxy/ks  across  each  pair  of  parallel  sides,  so  the  corresponding  shear 
strain  averaged  over  the  block  is  2zxy/wks  .  The  net  shear  compliance  is  thus 

1  1  ^  2 

G12  G  wks 

The  joint  normal  stiffness  is  the  only  original  elastic  constant  that  varies  with  stress.  For  this 
analysis  we  will  hold  it  constant  at  its  value  when  joint  closure  is  half  of  joint  thickness. 


Next,  we  must  solve  the  characteristic  equation 
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which  is  a  necessary  condition  for  the  general  functional  form  of  the  complex  potentials  assumed 
by  Lekhnitski  to  satisfy  the  equation  of  strain  compatibility  expressed  in  terms  of  stresses.  The 
two  positive  roots  can  be  written 
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where,  after  using  the  above  relations  among  elastic  constants, 

2  1  1  -  v  +  2 G/wkN 

2  1  -  v  +  2G/wks 


The  essence  of  Lekhnitski’s  procedure  is  that  if  $>j(z) ,  <J>2(z)  are  practically  any  complex 
valued  potential  functions  of  a  complex  variable  z,  then  the  stress  fields 
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ox-  -2Re^(Zj)  +  ^2(22)] 

ay  -  2Re[0/1(21)+a>2(22)]  ft 

V  -  2Im[X10/1(21)  +  X20/2(z2)] 

will  satisfy  all  the  equations  of  two-dimensional  static  orthotropic  elasticity  provided 

zk  -  x  +  (Kky  ,  £  -  1,2;  i  -  (-1)1/2  •  * 


The  task  now  is  to  find  potential  functions  for  which  the  stresses  satisfy  the  appropriate 
boundary  conditions.  For  remote  uniaxial  stress  ah  in  the  horizontal  (x)  direction  the 
appropriate  functions  are 

/  flo.  zk  +  [z l  +  a  2(X.  -  1)]1/2 

-  — -r  -  ,  *  -  .  -WAA -  '  *-‘'2 

where  a  is  the  radius  of  the  opening.  This  completes  the  derivation  of  the  fundamental  solution. 


The  problem  of  interest  is  that  of  a  tunnel  loaded  by  vertical  and  horizontal  free-field 
stresses  as  well  as  internal  loads  due  to  the  liner.  We  will  approximate  the  latter  by  a  uniform 
internal  pressure  p{  .  The  full  solution  for  any  stress  component  ai  can  be  built  up  by 
superposition  as  follows.  Let  d^r.to)  be  the  stress  component  at  radius  r  from  the  center  of  the 
tunnel  and  angle  to  from  the  horizontal,  corresponding  to  a  unit  remote  load  in  the  horizontal 
direction.  Then  by  linearity,  the  solution  for  horizontal  free  field  stress  oh  will 
be  ah  dfato)  .  This  corresponds  exactly  to  the  fundamental  solution  derived  above.  For  the 
problem  of  interest  there  are  two  other  sources  of  loading,  vertical  free  field  stress  ov  and 
internal  pressure  p{.  Because  the  implicitly  jointed  medium  is  symmetric  about  45°  lines,  the 
solution  for  the  former  is  ov&£r,ji/ 2 -to).  For  the  internal  pressure  we  superpose  three  more 
fields:  a  uniform  pressure  pt  everywhere,  and  the  fields  6((r,co)  ,  -p/6((r,ji/2-a>)  due  to 
horizontal  and  vertical  free  field  tensions  to  cancel  out  the  uniform  pressure  there.  This 
completes  the  solution  for  the  stresses  around  a  statically,  biaxially  loaded,  internally  pressurized 
opening  in  an  infinite  orthotropic  elastic  medium. 


To  apply  this  analysis  to  the  situation  at  peak  stress  arrival  time,  we  set 
ov  -  100  MPa,  ah  «•  33  MPa,  and  pt  -  8  MPa.  One  final  adjustment  is  needed:  because 
as  already  noted  the  free  field  stress  varies  in  the  vertical  direction,  after  evaluating  the  stresses 
by  superposition  as  described  above,  the  result  is  then  scaled  down  by  an  amount  representing 
the  locally  lower  free  field  stress  level.  The  results  are  plotted  in  Figures  4-13  and  4-14  as  the 
curves  labelled  "elastic".  The  agreement  with  numerically  derived  results  is  striking. 
Disagreements  arise  where  expected,  i.e.,  when  stresses  are  high  and  plasticity  comes  strongly 
into  play  (radiai  stresses  at  ±45°  at  intermediate  radii,  tangential  at  0°  and  ±90°  at  the  tunnel 
wall).  But  the  intermediate  maxima  in  tangential  stress  along  the  ±45°  radials  are  mimicked,  and 
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decay  rates  beyond  the  maxima  in  all  cases  are  well  approximated.  The  radial  stres  js  above  and 
below  the  tunnel  (Figure  4-13)  closely  agree  with  the  numerical  solutions.  This  confirms  the 
earlier  statements  concerning  the  extent  of  the  relief  provided  by  the  tunnel  wall. 

In  Figure  4-14  the  vertical  scale  for  tangential  stress  was  cut  off  at  160  MPa  for  clarity, 
owing  to  the  very  large  stress  concentrations  in  the  elastic  solution.  In  particular,  the  tangential 
stress  at  the  tunnel  wall  along  the  horizontal  radial  comes  out  at  673  MPa.  This  compares  with 
a  value  of  259  MPa  from  an  isotropic  elastic  analysis  under  the  same  loading  conditions.  In  the 
orthotropic  elastic  solution  the  high  stress  concentration  is  very  localized  near  the  tunnel  wall. 
In  the  first  meter,  the  tangential  stress  drops  to  144  MPa,  i.e.,  by  a  factor  of  4.7.  As  such,  the 
orthotropic  elastic  solution  would  not  be  expected  to  accurately  represent  the  stresses  near  the 
tunnel  in  an  explicitly  jointed  medium,  even  if  the  blocks  and  joints  remained  elastic.  However, 
the  more  severe  concentration  under  homogeneous  orthotropic  conditions  (compared  with 
isotropic)  suggests  that  the  global  orthotropy  of  the  jointed  medium  may  be  leading  to  more 
severe  tunnel  deformations  than  would  arise  in  an  isotropic  medium  with  the  same  vertical  and 
horizontal  properties. 

4.3.4  Remarks  on  the  Numerical  Solutions  to  Problem  4. 

Three  of  the  five  participants  (CRT,  WA,  Itasca)  obtained  numerical  solutions 
which — where  comparisons  were  possible — agree  with  each  other  in  most  practical  respects.  All 
of  these  solutions  appear  credible,  based  on  the  significant  body  of  evidence  available,  i.e., 

•  Use  of  rational  continuum  models  to  represent  the  rock, 

•  Use  of  physically  based  models  for  the  joints, 

•  Compatibility  of  results  with  basic  understanding  of  wave  propagation  processes, 

•  Absence  of  obvious  numerical  artifacts  such  as  spurious  reflections, 

•  Comparison  of  stresses  and  strains  with  complete  and  partial  analytic  solutions  in 
all  the  benchmark  problems. 

RE/SPEC’s  numerical  solution  to  problem  4  differs  significantly  from  the  first  three,  for 
example  in  springline  closure  (Figure  4-12),  radial  stresses  above  and  below  the  tunnel  at  peak 
stress  arrival  (Figure  4-13),  and  radial  and  tangential  stresses  along  the  other  radials  at  the  end 
of  the  positive  phase  (Figures  4-15  and  4-16).  In  the  case  of  the  stresses  at  the  peak  stress  arrival 
time,  RE/SPEC’s  solution  also  differs  from  the  orthotropic  elastic  solution,  while  the  latter  agrees 
well  with  CRT  and  Itasca’s  numerical  results  and  does  so  for  reasons  that  we  believe  we 
understand.  The  large  springline  closure  could  be  related  to  a  spurious,  compressive  reflection 
from  the  lower  boundary  of  the  computational  grid.  There  is  no  obvious  explanation  for  the 
significant  residual  stresses  around  the  tunnel  at  the  end  of  the  positive  phase  (Figures  4-15  and 
4-16).  There  do  appear  to  be  explanations  for  most  of  the  other  identified  shortcomings  and 
breakdowns  in  RE/SPEC’s  numerical  solutions  to  the  earlier  problems,  and  with  further  effort  the 


results  could  probably  be  corrected.  However,  the  totality  of  evidence  developed  during  this 
study  suggests  that  this  solution  to  problem  4  is  not  as  accurate  as  the  aforementioned  three. 

LLNL’s  approach  omits  a  great  deal  of  pertinent  physics,  most  glaringly  the  behavior  of 
the  rock  as  a  deformable  continuum.  A  glance  at  Figure  4-19  shows  convincingly  that  when  the 
blocks  as  specified  in  problem  4  are  permitted  to  deform,  they  will  do  so,  and  their  deformations 
will  contribute  substantially  to  the  conventional  indices  of  tunnel  distress,  viz.,  tunnel  closure. 
Moreover,  particularly  with  respect  to  springline  closure,  the  coupling  between  vertical  and 
horizontal  stresses  and  deformations  plays  a  crucial  role.  There  is  no  mechanism  in  the  LLNL 
model  for  directly  representing  this  coupling.  LLNL  attempted  to  approximate  the  effect  by 
imposing  a  constant  33  MPa  horizontal  stress.  The  result  is  seen  in  Figure  4-12:  even  with  the 
imposed  horizontal  stress,  the  springling  closure  is  still  negative  (i.e.,  the  springline  elongates), 
in  contrast  with  the  expected  result  in  the  face  of  dilatant  rock  response.  More  generally,  the 
pattern  of  deformations  depends  crucially  on  the  distribution  of  stresses  near  the  tunnel,  and  the 
LLNL  approach  cannot  accurately  estimate  these  stresses.  It  is  for  these  reasons  that  the 
approach  seems  ill-suited  for  the  problems  at  hand.  This  is  not  to  say  that  there  aren’t  certain 
problems  for  which  it  may  be  applicable.  But  to  be  useful  it  would  have  to  be  accompanied  by 
guidelines  for  its  applicability. 


SECTION  5 


CONCLUSIONS 


The  primary  objective  of  the  benchmark  calculation  study  was  to  understand  the  effect 
of  computational  approach  on  predictions  of  tunnel  deformations.  We  can  only  claim  success 
in  meeting  that  objective  by  adopting  a  rather  broad  view  of  what  comprises  a  computational 
approach.  That  is,  we  must  include  more  than  purely  numerical  issues.  In  most  cases  where  we 
have  identified  problems  or  differences  and  where  there  is  an  apparent  cause,  that  cause  lies  in 
the  way  the  medium  was  idealized.  For  example, 

(1)  The  LLNL  approach  oversimplified  the  continuum  physics  of  the  rock  by  lumping 
it  in  the  joints,  and  this  probably  led  to  their  substantially  different  results  from 
all  the  others. 

(2)  In  some  early  solutions  by  WA  the  normal  stiffness  of  the  joints  was  based  on  the 
properties  of  the  intact  rock  and  therefore  unrepresentative  of  the  specified  joint 
stiffness.  The  code  was  modified  so  that  a  constant  stiffness  approximating  the 
actual  normal  joint  stiffness  could  be  specified.  The  joints  still  did  not  have  the 
specified  nonlinear  behavior,  but  they  did  have  a  constant  stiffness  much  closer 
to  the  average  over  the  expected  range  of  joint  closures.  With  this  modification, 
the  reason  for  the  stiffness  deviation  was  understood,  and  was  due  to  a  conscious 
choice  in  an  area  that  could  now  be  regarded  as  material  idealization. 

(3)  Itasca  and  WA  both  chose  to  use  an  isotropic  model  for  the  implicitly  jointed 
region.  A  possible  result  of  this  is  the  discrepancy  in  horizontal  stresses  across 
the  implicit/explicit  boundary  noted  in  Figure  4-9. 

There  are  also  a  few  examples  of  issues  more  clearly  and  classically  numerical  in  nature: 

(4)  The  WA  sliding  interface  model  was  causing  a  great  deal  of  numerical  chatter  in 
problem  2,  the  first  one  in  this  study  where  it  was  used.  After  changing  certain 
features  of  the  model  other  than  the  stiffness  (see  Section  E4),  the  problem  was 
resolved. 

(5)  RE/SPEC  had  a  problem  with  the  nonreflecting  boundary  in  problem  3,  and 
recognized  that  it  was  caused  by  the  wave  speeds  used  in  the  boundary  condition. 
They  did  not  elect  to  remedy  the  problem  during  the  course  of  this  study.  The 
spurious  bottom  reflection  probably  caused  higher  stresses  at  the  end  of  the  free 
field  positive  phase  (compared  with  other  participants),  and  this  in  turn  may  have 
led  to  significantly  larger  springline  closures. 

It  is  important  to  recognize  that  there  is  no  definitive  proof  of  the  cause-and-effect 
relationships  suggested  above.  They  rely  more  on  physical  reasoning  and  induction  than 


deduction.  In  fact,  as  indicated  for  example  in  the  discussion  on  divergent  flow,  glib  physical 
arguments  may  either  fail  or  reveal  themselves  as  oversimplifications  when  placed  under  closer 
scrutiny. 

In  particular,  statement  (3)  above,  while  apparently  logical,  could  be  questioned.  As 
suggested  in  Section  4.2.3,  by  judicious  choice  of  material  parameters,  an  isotropic  elastic  model 
could  be  made  to  give  the  same  coupling  between  vertical  and  horizontal  stresses  and  strains  as 
an  orthotropic  one.  The  material  parameters  quoted  in  Appendix  E  for  WA’s  implicit  model 
correspond  to  a  smaller  Poisson’s  ratio  than  for  intact  rock,  as  do  the  anisotropic  elastic 
relationships  of  Section  4.3.3,  when  used  with  the  unstrained  stiffness  of  the  joints.  But  as  the 
joints  compress  and  stiffen,  Poisson’s  ratio  increases.  The  situation  is  obviously  more  complex 
than  first  meets  the  eye. 

As  to  the  overall  success  of  the  calculations  themselves,  the  remarks  in  Section  4.3.4 
summarize  the  author’s  evaluation.  But  even  the  organizations  which  appeared  to  successfully 
negotiate  problem  4  did  not  do  so  without  breakdowns  and  detours  along  the  way.  If  the  tunnel 
deformation  calculations  had  not  been  preceded  by  a  sequence  of  progressively  more  challenging 
problems,  the  results  would  have  been  much  more  scattered.  Often  a  calculator  is  faced  with  a 
more  complex  problem  than  any  of  those  treated  here,  and  fewer  resources  to  devote  to  it.  Some 
of  the  lessons  are  obvious:  the  calculator  should  take  care  to  include  all  the  correct  physics  he 
can,  while  the  customer  should  support  as  much  preliminary  work  as  possible,  and  barring  that, 
maintain  a  healthy  sense  of  skepticism  about  the  outcome  of  any  isolated,  complex,  numerical 
simulation. 
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APPENDIX  A 


DESCRIPTION  OF  THE  PR0NT02D  COMPUTER  CODE  USED  IN  THE 
UNDERGROUND  TECHNOLOGY  PROGRAM  BENCHMARK  ACTIVITY 

John  Osnes  (RE/SPEC),  Ben  H.  Thacker  (SwRI), 
and  David  S.  Riha  (SwRI) 


Al.  INTRODUCTION. 

PR0NT02D  is  a  two-dimensional  transient  solid  dynamics  code  for  analyzing  large 
deformations  of  highly  nonlinear  materials  subjected  to  extremely  high  strain  rates  (Taylor,  L.  M. 
and  Flanagan,  D.  P.,  1987).  It  is  the  latest  in  a  series  of  transient  dynamics  finite  element  codes 
that  has  been  developed  at  Sandia  National  Laboratories,  beginning  with  HONDO  (Key,  S.  W., 
et  al.,  1978).  As  such,  PR0NT02D  contains  a  number  of  state-of-the-art  numerical  algorithms, 
including  an  adaptive  timestep  control  algorithm,  a  robust  hourglass  control  algorithm,  a  very 
accurate  incremental  rotation  algorithm,  and  a  robust  surface  contact  algorithm.  Four-noded, 
uniform-strain,  quadrilateral  elements  with  single-point  integration  are  used  in  the  finite  element 
formulation.  Beyond  its  general  capabilities,  PR0NT02D  was  chosen  for  the  benchmark 
calculations  because  new  constitutive  models  are  readily  added  and  because  a  three-dimensional 
derivative  of  the  program  is  available  (PR0NT03D)  if  needed  by  the  UTP  in  the  future. 

The  two  features  that  make  the  UTP  benchmark  calculations  relatively  unique  among 
large-deformation  solid  dynamics  problems  are: 

(1)  A  jointed  rock  mass  subject  to  loading  conditions  that  could  result  in  large  relative 
motions  (sliding)  between  adjacent  rock  blocks,  and 

(2)  Loading  conditions  that  could  result  in  substantial  tensile  failure,  plastic 
deformation,  and  dilation  within  the  rock  blocks. 

The  joints  also  affect  the  elastic  behavior  of  the  rock  mass  by  reducing  its  elastic  moduli  from 
the  intact  rock  values.  Further,  the  elastic  moduli  of  jointed  rock  generally  are  nonlinear 
functions  of  the  joint  apertures. 

The  surface  contact  algorithm  in  PR0NT02D  is  directly  applicable  to  modeling  the  slip 
between  adjacent  rock  blocks.  It  is  designed  to  simulate  smooth,  cohesionless  surfaces  with  shear 
strengths  defined  in  terms  of  static  and  dynamic  coefficients  of  friction.  In  the  benchmark 
problems,  the  static  coefficient  of  friction  is  the  tangent  of  the  joint  friction  angle.  The  friction 
coefficient  does  not  decay  during  sliding,  so  the  dynamic  coefficient  of  friction  is  equal  to  the 
static  coefficient. 

The  surface  contact  algorithm  does  not  provide  a  means  to  model  the  elastic  behavior  of 
the  joints.  Consequently,  the  approach  used  in  the  benchmark  problems  is  to  superpose  the 
nonlinear  stress-displacement  response  of  the  joints  and  the  linear  stress-strain  response  of  the 
intact  rock.  The  result  is  a  composite  material  whose  elastic  behavior  is  equivalent  to  a  jointed 
rock  mass.  This  approach  has  been  developed  and  used  successfully  by  other  modelers 
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(Labreche,  D.  A.  and  Petney,  S.  V.,  1987),  and  it  is  referred  to  as  the  Complaint  Joint  Model 
(CJM).  Note  that  the  elastic  behavior  of  jointed  rock  masses  as  simulated  by  the  CJM  is 
inherently  anisotropic  because  the  stress-strain  response  normal  to  a  joint  set  is  substantially 
different  from  the  response  tangential  to  the  joint  set. 

Nine  constitutive  models  are  available  in  PR0NT02D.  They  include  an  elastic-plastic 
model  with  a  Mises  yield  criterion  that  is  directly  applicable  to  modeling  the  tunnel  liner.  The 
models  also  include  an  elastic-plastic  model  with  a  two-invariant  yield  surface.  However,  the 
flow  rule  in  that  model  is  nonassociative  and  cannot  be  reduced  to  an  associative  form. 
Consequently,  to  perform  the  benchmark  calculations,  a  new  constitutive  model  (Drucker-Prager) 
was  implemented  in  PR0NT02D  for  modeling  the  plastic  deformation  within  the  rock  blocks. 
The  new  model  also  simulates  the  limited  tensile  strength  of  the  intact  rock.  The  CJM  is 
incorporated  within  the  model  to  simulate  the  elastic  response  of  the  intact  rock,  including  joint 
sets. 


In  summary,  the  approach  to  modeling  the  jointed  rock  mass  uses: 

(1)  The  surface  contact  algorithm  in  PR0NT02D  to  explicitly  model  slip  between 

rock  blocks  along  joints. 

(2)  A  new  constitutive  model  that: 

(a)  Models  the  elastic  behavior  of  the  rock  and  joint  sets  as  an  equivalent 
composite  material  using  the  Compliant  Joint  Model  to  calculate  the 
nonlinear,  anisotropic  stress-strain  response  of  the  rock  mass. 

(b)  Simulates  yield  of  the  intact  rock  using  an  incremental  plasticity  method 
with  a  Drucker-Prager  yield  function  and  the  associated  plastic  potential 
function. 

(c)  Simulates  tensile  failure  using  a  limited-tension  algorithm  to  reduce 
principal  stress  components  that  exceed  the  tensile  strength  to  the  tensile 
strength  while  maintaining  the  original  principal  stress  directions. 

A2.  NUMERICAL  FORMULATION. 

A2.1  Explicit  Time  Integration. 

The  equations  of  motion  are  integrated  using  a  modified  central  difference  scheme  in 
PR0NT02D.  The  velocities  are  integrated  with  a  forward  difference  and  the  displacements  are 
integrated  with  a  backward  difference.  This  scheme  is  expressed  as 


fEXT  fINT 
Jt  ~  Jt 


“/♦*  ’  Ut  +  ^tat  +  At 


where  ftEXT  is  the  external  nodal  force,  ftINT  is  the  internal  noda1  force,  m  is  the  nodal  point 
lumped  mass,  and  At  is  the  time  increment.  This  central  differei.cc  operator  is  conditionally 
stable  and  the  Courant  stability  limit  is  given  by  the  highest  eigenvalue  of  the  system 
(Bathe,  K.  J.  and  Wilson,  E.  L.,  1976). 


Flanagan  and  Belytschko  (1984)  provided  eigenvalue  estimates  for  the  uniform  strain 
quadrilateral  used  in  PR0NT02D. 


Numerical  damping  is  introduced  in  the  solution  by  adding  artificial  viscosity.  This 
prevents  high  velocity  gradients  from  collapsing  an  element  before  is  has  a  chance  to  respond 
and  to  quiet  truncation  frequency  ringing.  The  technique  used  in  PRONT02D  is  to  add  viscosity 
to  the  "bulk"  response.  This  generates  a  viscous  pressure  in  terms  of  the  volume  strain  rate. 

A2.2  Four  Node  Uniform  Strain  Element. 


PR0NT02D  uses  a  four-noded  two-dimensional  uniform  strain  element  in  the  finite 
element  formulation.  A  one  point  integration  of  the  element  under-integrates  the  element  but 
provides  a  large  computational  advantage  over  a  two-by-two  integration  rule.  However,  this 
results  in  a  rank  deficiency  for  the  element  that  may  cause  spurious  zero  energy  (hourglass) 
modes. 


The  mean  stress-strain  formulation  for  the  uniform  strain  element  considers  only  a  fully 
linear  field.  Any  remaining  nodal  velocity  field  is  the  hourglass  field.  Possibly  severe 
unrestricted  mesh  distortion  can  occur  if  these  modes  are  excited.  The  method  used  in 
PR0NT02D  isolates  the  hourglass  modes  so  they  may  be  treated  independently  of  the  rigid  body 
and  uniform  strain  modes  (Flanagan,  D.  P.  and  Belytschko,  T.,  1981). 

A2.3  Material  Behavior. 

Several  classical  yield  surfaces  that  are  defined  in  terms  of  the  first  stress  invariant  and 
the  second  deviatoric  stress  invariant  have  been  used  to  model  plasticity  in  rock.  They  include 
the  linear  Mohr-Coulomb  and  the  Drucker-Prager  yield  criteria.  These  criteria  have  been 
reviewed  by  Callahan  and  Fossum  (1982).  The  Drucker-Prager  yield  function  was  used  for 
modeling  the  plasticity  in  the  benchmark  problems.  This  yield  function  is  defined  as 
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where  o-  are  the  components  of  the  stress  tensor,  a  and  k  are  material  constants,  Ij  is  the  first 
invariant  of  the  stress  tensor,  and  J2  is  the  second  invariant  of  the  deviatoric  stress  tensor.  The 
yield  surface  is  the  locus  of  stress  states  at  which  the  value  of  the  yield  function  is  zero  (f  =  0). 
Consequently,  the  Drucker-Prager  yield  surface  is  a  cone  in  principal  stress  space  with  the  axis 
of  the  cone  along  the  hydrostat.  The  specification  of  perfect  plasticity  in  the  benchmark 
problems  means  that  the  Drucker-Prager  yield  surface  does  not  change  as  a  result  of  plastic 
deformation  (i.e.,  the  values  of  a  and  k  do  not  change). 

The  Drucker-Prager  and  the  Mohr-Coulomb  yield  criteria  are  equivalent  only  at  certain 
stress  states,  depending  on  how  the  Drucker-Prager  material  constants  are  evaluated.  To  match 
the  Mohr-Coulomb  criterion  in  triaxial  compression,  the  material  constants  must  be  evaluated  as 
follows: 


a  - 


2  sin 


v/3~(3  -sin<j>) 


k  - 


2^5"  COS(j> 


y  3  -  sin  0  ) 


where  <(>  and  cQ  are  the  Mohr-Coulomb  friction  angle  and  cohesion,  respectively.  The  resultant 
Drucker-Prager  yield  surface  circumscribes  the  Mohr-Coulomb  yield  surface. 


The  stress  state  is  elastic  when  the  yield  function  is  negative  (f  <  0).  When  the  stress 
state  is  on  the  yield  surface  (f  -  0)  and  the  loading  condition  is  such  that  it  would  result  in  fit  0, 
the  resulting  deformation  is  plastic.  The  components  of  the  resultant  plastic  strain  tensor  are 
defined  in  terms  of  the  flow  rule,  which  is  classically  stated  as  proportional  to  th„  gradient  of  a 
plastic  potential  function  g(0;j).  The  proportionality  constant  is  determined  by  the  consistency 
condition  which  requires  that  the  stress  state  remain  on  the  yield  surface  during  plastic 
deformation.  When  the  plastic  potential  function  and  the  yield  function  coincide  for  all  stress 
.’tates  (g  =  f),  the  flow  rule  is  said  to  be  associated  with  the  yield  function,  as  required  in  the 
oenchmark  problems. 

An  incremental  method  using  tangent  stiffness  (Chen,  W.  F.  and  Han,  D.  J.,  1988)  is  used 
to  implement  Drucker-Prager  plasticity  in  the  cor  ttitutive  model  that  was  added  to  PR0NT02D 
for  the  benchmark  problems.  The  incremental  stress-strain  relationship  can  be  expressed  in  the 
following  general  form: 
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where  da ^  and  dzkl  are  the  components  of  the  incremental  stress  and  strain  tensors,  respectively, 
Cijkl  are  the  components  of  the  elastic  coefficient  tensor,  and  repeated  subscripts  indicate 
summation.  The  coefficient  tensor  in  the  brackets  represents  the  elastic-plastic  tensor  of  tangent 
moduli  for  an  elastic-perfectly  plastic  material.  The  quantities  in  the  brackets,  including  C-kl  for 
a  nonlinear  elastic  material  such  as  in  the  CJM,  are  evaluated  at  the  current  stress  state. 

For  the  Drucker-Prager  yield  surface,  the  partial  derivatives  of  the  yield  function  with 
respect  to  the  stress  components  are 


d 


where  s-  are  the  components  of  the  deviatoric  stress  tensor  (s-  =  o- -  6^/ j/3)  and  6,y  is  the 
Kronecker  delta  function.  In  the  constitutive  model  added  to  PR0NT02D  the  following  extended 
form  of  the  plastic  potential  function  associated  with  the  Drucker-Prager  yield  function  is  used: 
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and  ip  is  the  dilation  angle.  The  plastic  potential  function  results  in  an  associated  flow  rule  when 
the  dilation  angle  equals  the  friction  angle  (so  that  B  =  a  and  g  -  f).  The  components  of  the 
flow  rule  (the  partial  derivatives  of  g  with  respect  to  o-)  are 
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A2.4  Limited-tension  Algorithm. 

A  limited-tension  algorithm  developed  by  Callahan  (Callahan,  G.  D.  and  Fossum,  A.  F., 
1982)  and  implemented  in  SPECTROM-32  (Callahan,  G.  D.,  et  al.,  1989)  is  used  in  the 
constitutive  model  to  simulate  tensile  failure.  In  this  algorithm  the  tensile  strength  at  each 
integration  point  Tx  is  set  initially  to  the  intact  tensile  strength  To-  At  each  integration  point  at 
each  time  step,  the  principal  stresses  (a}  au  *  ojn)  and  their  unit  direction  vectors 
(  rij,  A jj,  nJU  )  are  calculated  based  on  the  trial  stress  state  o »  at  that  location  and  time.  Each 
principal  stress  component  that  exceeds  the  tensile  strength  at  the  integration  point  Tx  is  set  to 
the  residual  tensile  strength  Tr .  In  addition,  the  tensile  strength  at  the  integration  point  is  set  to 
Tr  if  Tx  is  exceeded.  (Recall  that  tension  is  negative  according  to  the  sign  convention,  so 
exceeding  Tx  implies  that  the  stress  is  less  than  -Tx  and  the  stress  is  set  to  -Tr  .)  The  resultant 
stress  state  is  calculated  by  transforming  the  principal  stress  components  back  to  the  global 
coordinate  system  while  maintaining  the  original  principal  stress  directions.  This  transformation 
is  accomplished  according  to  the  following  equation: 


°ij  "  niinij°i  +  nmnujaii  +  nuiinnijaiii 

where  nfi  is  the  Ith  component  c c  unit  direction  vector  h} . 

In  the  limited-tension  algorithm,  the  stress  state  is  modified  in  a  very  direct  way  when 
the  tensile  strength  is  exceeded.  The  deformation  field  also  is  affected  because  modifying  the 
stress  state  produces  an  imbalance  between  the  external  forces  and  the  internal  forces  (the  integral 
of  the  stress  divergence).  In  turn,  this  imbalance  results  in  accelerations  and  ultimately 
deformation  in  response  to  the  tensile  failure.  Changing  the  tensile  strength  to  a  residual  value 
after  tensile  failure  allows  the  formation  of  tension  cracks  to  be  simulated  by  a  reduction  in 
strength.  However,  this  approach  to  simulation  of  tensile  cracking  is  very  crude  because  it  does 
not  account  for  crack  orientation  and  the  resultant  anisotropy  in  the  tensile  strength.  It  is 
justifiable  when  the  principal  stress  directions  are  fairly  constant  throughout  a  simulation,  so  that 
the  tensile  direction  and  the  resultant  orientation  of  tension  cracks  are  fairly  constant.  In  general, 
this  condition  could  not  be  ensured  in  the  benchmark  calculations.  Consequently,  the  tensile 
strength  is  maintained  at  the  intact  value  even  after  tensile  failure  (i.e.,  Tr  =  7^). 

A2.5  Implicit  Rock  Joint  Modeling. 

The  compliant  joint  model  (CJM)  is  a  three-dimensional  elastic  model  used  to  represent 
the  mechanical  response  of  a  rock  matrix  with  joints  (Labreche,  D.  A.  and  Petney,  S.  V.,  1987). 
The  model  implemented  in  PR0NT02D  is  a  two-dimensional  version  of  the  CJM.  The  CJM 
assumes  that  the  deformation  of  a  jointed  rock  mass  is  caused  by  the  combination  of  the  elastic 
response  of  the  rock  matrix  and  the  elastic  response  of  the  joints.  The  unfractured  rock  matrix 
is  modeled  as  an  isotropic  linear  elastic  solid  and  the  joints  as  a  nonlinear  elastic  layer.  The 
model  consists  of  up  to  four  joint  sets  with  various  spacing  and  orientation.  The  response  of  the 
CJM  is  effectively  anisotropic  due  to  the  presence  of  the  joint  sets. 

The  displacement  behavior  normal  to  the  joint  plane  varies  nonlinearly  with  stress  and  the 
response  is  governed  by  the  physical  characteristics  such  as  joint  aperture,  roughness  and  contact 
between  joint  surfaces,  strength  of  the  wall  rock,  and  presence  or  absence  of  filling.  The 
response  is  also  governed  by  the  recent  movement  history  of  the  rock  joint.  The  aperture  change 
of  a  well-mated  joint  without  filling  can  be  represented  by  a  hyperbolic  relationship  between  the 
stress  normal  to  joint  plane  and  the  joint  closure  (Bandis,  S.  C,  et  al.,  1983  and  Goodman,  R.  E., 
1976).  The  form  used  in  the  CJM  is  given  by  the  equation 


where, 

an  -  normal  stress  acting  on  the  joint  (compression  positive) 
a  -  half-closure  stress  (stress  to  reduce  the  aperture  to  6/2 
vw  -  joint  normal  displacement  (closure  positive,  zero  when  o;j  =  0) 
6  -  maximum  joint  closure 

m  -  exponent 


The  stiffness  normal  to  the  joint,  kn  ,  is  defined  by 
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which  is  the  slope  of  the  nonlinear  joint  closure.  For  the  Benchmark  problems,  a  -  125  MPa  and 
m-  1.  This  hyperbolic  function  is  used  in  the  numerical  formulation  for  the  normal 
stress-displacement  response  when  the  normal  stress  is  compressive.  A  constant  value  of  10*4 
times  the  dilatational  modulus,  X  2^,  of  the  matrix  is  used  for  joint  normal  stiffness  when  the 
normal  stress  is  tensile.  This  stiffness  defines  a  linear,  normal  stress-normal  displacement 
telationship  for  the  joint  in  tensior .  This  small  but  nonzero  stiffness  avoids  numerical 
singularities  in  the  computations.  The  transmission  of  small  tensile  stresses  by  elements  having 
moduli  four  orders  of  magnitude  lower  than  the  matrix  are  accepted  as  insignificant  relative  to 
typical  compressive  stresses. 

The  compressive  normal  stress-displacement  relation  is  nonlinear  and  elastic.  Thus,  there 
is  no  path  dependence  and  all  joint  closure  is  fully  recoverable  in  this  model. 

The  shear  response  of  joints  generally  consists  of  a  nearly  linear  rise  to  peak  shear 
strength,  followed  by  softening  behavior  with  continuing  shear.  However,  it  is  acknowledged  that 
joint  size  dependence  may  alter  this.  The  peak  shear  strength  response  of  joints  is  supported  by 
studies  such  as  Goodman  (1976)  and  Bandis  (Bandis,  S.  C.,  et  al.,  1983).  The  shear  stress  of 
joints  is  known  to  be  inelastic  when  the  stress  state  approaches  the  peak  shear  strength.  The 
CJM  only  models  the  elastic  joint  shear  response.  The  shear  stiffness,  which  is  independent  of 
the  normal  stiffness  across  ihe  joint,  relates  the  shear  displacement  to  shear  stress. 


where, 

x  -  shear  strength 

ks  -  joint  shear  stiffness 

v(  -  joint  displacement  tangent  to  plane  of  joint 


The  shear  stiffness  is  defined  by 


dx 


dvt 


The  CJM  introduces  anisotropy  by  reducing  the  composite  modulus  of  the  modeled  rock 
mass  in  the  direction  normal  to  the  joint  in  an  otherwise  linear  elastic  material.  Computation  of 
effective  moduli  is  used  to  implement  the  softening  effects  of  the  joints.  Stress  equilibrium  is 
the  underlying  assumption  in  the  combined  response  of  the  rock  matrix  and  the  joint  sets.  The 
normal  stress  across  the  joint  is  equal  to  the  stress  in  the  rock  matrix  normal  to  the  joint  plane. 
The  shear  stress  on  the  joint  is  equal  to  the  shear  stress  in  the  rock  matrix  parallel  to  the  joint 
plane. 


Correspondingly,  the  strains  in  the  rock  matrix  are  combined  with  the  strains  in  the  joints 
additively.  That  is,  the  effective  compliance  for  the  rock  mass  is  the  sum  of  the  compliances  of 
the  unfractured  rock  and  the  joint  sets.  The  effective  Young’s  modulus,  E  ,  normal  to  a  joint 
and  the  effective  shear  modulus,  G  ,  parallel  to  a  joint,  are 

1  _  1  ^  1 

E  E  kns 

J.  -  JL  +  _L 

G  ~G  T^s 

where  E  and  G  are  the  elastic  moduli  of  the  rock  matrix,  kn  and  ks  are  the  normal  and  shear 
stiffness  of  the  joint  set,  and  s  is  the  spacing  of  the  joints  in  the  joint  set. 

The  description  of  the  compliant  joint  model  is  for  an  elastic  material  that  models  the 
elastic  response  of  a  jointed  rock  mass.  However,  the  implementation  of  the  limited  tensile 
algorithm  and  Drucker-Prager  plasticity  allows  a  more  realistic  representation  of  the  rock  matrix. 

A2.6  Explicit  Rock  Joint  Modeling. 

PR0NT02D  treats  contact  as  a  kinematic  constraint.  That  is,  the  algorithm  modifies  the 
accelerations  of  the  nodes  along  the  contact  region  such  that  the  kinematic  constraints  are 
satisfied.  The  algorithm  uses  a  partitioned  approach  to  enforce  compliance  between  two  contact 
surfaces  by  allowing  each  surface  to  act  as  a  master  surface  for  a  fraction  of  each  time  step  and 
a  slave  for  the  remainder. 

There  are  four  steps  involved  in  the  contact  algorithm.  First,  the  contact  surface  geometry 
is  recalculated  at  each  time  step  and  the  predicted  configuration  is  computed  by  integrating  the 
motion  without  regard  to  the  kinematic  constraints  required  by  the  contact  surfaces.  The 
following  quantities  are  calculated  for  each  node: 

a  «  flm 
v  «  v  +  A  td 
x  -  x  +  A  tv 

where  /  is  the  residual  force  vector  (sum  of  external  forces  minus  the  sum  of  internal  forces),  m 
is  the  nodal  mass,  v,  is  the  current  velocity,  and  At  is  the  time  increment.  The  predicted 
kinematic  quantities  are  denoted  by  the  hat 

The  second  step  is  surface  tracking,  or  the  process  of  matching  nodes  along  one  surface 
with  the  mating  surface.  The  algorithm  used  is  to  locate  the  spatially  nearest  master  node  to  the 
possible  point  of  contact.  This  procedure  can  be  the  most  time  consuming  portion  of  the  analysis 
for  this  type  of  problem.  To  streamline  the  tracking  algorithm,  the  nearest  master  node  to  a 
given  slave  node  at  one  time  step  is  assumed  to  be  in  the  vicinity  of  the  nearest  master  node  at 
the  next  time  step.  Therefore,  at  each  time  step,  the  nearest  master  node  is  the  starting  point  for 
the  search  along  the  surface. 


The  next  step  is  to  determine  contact  or  penetration.  The  slave  node  is  oriented  with 
respect  to  the  master  segments  connected  to  the  tracked  master  node.  The  depth  and  position 
coordinates  are  calculated  for  both  master  segments  connected  to  the  nearest  master  node.  From 
these  quantities,  if  the  depth  is  positive,  the  slave  node  is  penetrating  the  segment  and  if  the 
position  is  positive  then  the  slave  node  is  along  the  segment.  When  the  master  surface  forms  an 
outside  comer,  there  may  be  penetration  of  both  segments.  In  this  case,  the  algorithm  determines 
with  which  master  segment  the  slave  surface  is  more  strongly  in  contact.  One  limitation  of  the 
algorithm  is  that  it  can  not  detect  a  contact  surface  contacting  itself. 

The  final  step  in  the  contact  algorithm  is  to  calculate  the  penetration  forces  imposed  on 
the  master  surface  by  the  slave  surface  to  restore  kinematic  compliance.  These  forces  are 
calculated  as  a  fraction  of  the  forces  that  would  be  imposed  by  the  slave  nodes  if  the  master 
surface  was  rigid.  This  fraction  is  based  on  the  fraction  of  each  time  step  for  which  the  surfaces 
act  as  master  and  slave.  The  roles  are  reversed  for  the  remaining  fraction  of  the  time  step.  The 
accelerations  are  calculated  to  predict  the  response  of  the  master  surface  to  these  penetration 
forces  such  that  the  response  of  each  contacting  slave  node  is  constrained  by  its  master  nodes. 
The  principle  of  virtual  work  is  employed  to  define  the  accelerations  of  the  master  nodes  in 
response  to  the  penetration  forces.  When  friction  is  present,  the  relative  tangential  motion  of  the 
contacting  slave  nodes  is  resisted.  The  magnitude  of  the  tangential  force  exerted  on  the  master 
surface  on  a  slave  node  cannot  exceed  the  friction  force.  Thus  friction  adds  a  tangential 
acceleration  to  the  nodes  in  contact 

A2.7  Transmitting  Boundary  Conditions. 

A  transmitting  boundary  is  used  to  simulate  a  semi-infinite  domain  outside  the  boundary, 
where  the  wave  speeds  of  the  material  on  both  sides  of  the  boundary  are  the  same.  The  region 
exterior  to  the  boundary  is  replaced  with  an  energy-absorbing  boundary  condition  that  behaves 
as  if  the  energy  is  transmitted  across  the  boundary.  Thus,  no  energy  is  reflected  back  into  the 
interior  region. 

The  nonreflecting  boundary  is  implemented  in  PR0NT02D  according  to  a  technique 
proposed  by  Lysmer  (Lysmer,  J.  and  Kuhlemeyer,  R.  L.,  1969).  The  basic  idea  is  to  apply 
tractions  at  the  boundary  that  will  exactly  cancel  the  stresses  that  would  be  reflected  from  a  free 
surface.  Hence,  the  numerical  technique  involves  calculating  and  applying  the  following  tractions 
to  the  surface  of  the  nonreflecting  boundary 

-  -p 

-  -PVsat 

where  on  and  xs  are  the  normal  and  shear  tractions;  p,  V  ,  and  Vs  are  the  current  density, 
compressional  wave  speed,  and  shear  wave  speed  of  the  material  along  the  boundary; 
and  un  and  u(  are  the  normal  and  tangential  components  of  the  current  velocity  at  the 
boundary.  At  each  time  step,  PR0NT02D  updates  the  tractions  at  the  boundary  using  the  current 
density  and  effective  dilatational  and  shear  moduli  in  each  element  along  the  boundary. 

In  calculating  the  current  wave  speeds,  PR0NT02D  assumes  that  the  wave  speeds  are 
isotropic  and  independent  of  the  orientation  of  the  boundary.  Obviously,  the  joint  sets  in  a  rock 


mass  yield  an  anisotropic  medium  in  which  the  wave  speeds  depend  on  the  incident  direction. 
However,  the  anisotropy  inherent  in  the  CJM  is  not  accounted  for  in  the  wave  speed  calcination 
in  PR0NT02D.  Consequently,  the  wave  speeds  used  to  calculate  the  tractions  along  the 
nonreflecting  boundary  are  somewhat  in  error,  the  magnitude  of  v/hich  depends  on  the  orientation 
of  the  boundary  with  respect  to  the  principal  directions  of  the  anisotropy.  This  error  results  in 
a  partial  reflection  of  the  incident  stress. 
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APPENDIX  B 


THEORETICAL  ASPECTS  OF  THE  LLNL  DIBS 
DISCRETE  ELEMENT  CODE  USED  IN  THE  DNA/UTP 
BENCHMARK  EXERCISE 

Francois  Heuze 


Bl.  INTRODUCTION. 

Calculations  for  this  study  were  performed  with  the  LLNL  Discrete  Interacting  Block 
System  (DIBS)  numerical  model  (Walton,  O.  R.,  1982  and  Walton,  O.  R.,  1980)  using  a  sew 
"rounded  comer"  version  of  that  model.  DIBS  was  started  from  the  algorithms  of  earlier 
"distinct-element"  models  (Cundall,  P.  A.,  1974),  and  evolved  at  LLNL.  It  calculates  the  motion 
of  each  discrete  polygonal  block  as  it  responds  to  contact  forces,  boundary  and  applied  loads,  and 
gravity.  Arbitrary  (convex)  polygonal  shapes  are  allowed  for  each  particle  or  block.  Several 
simplifying  assumptions  make  the  calculations  tractable.  Ail  deformations  take  place  at  the 
contacts,  which  include  the  compliance  of  the  rock  blocks  as  well  as  that  of  the  rock  joints; 
meanwhile  the  blocks  are  kept  rigid  during  the  calculations.  Note  that  all  rock  joints  in  DIBS 
are  explicit;  there  is  no  implicit  jointing.  All  starting  contacts  between  polygons  are  assumed  to 
be  comer-on-side  (i.e.,  temporary,  sliding  joint-elements  or  contacts  are  set-up  at  polygonal 
vertices,  as  needed).  Small  but  finite  "overlaps"  occur  as  normal  forces  are  developed.  Similarly 
a  finite,  partially  recoverable  shear-strain  develops  in  the  contact  before  frictional  sliding  is 
initiated.  In  the  rounded-comer  version  of  the  model,  comer-comer  contacts  are  also  allowed, 
and  they  serve  to  allow  a  smooth  transition  in  the  direction  of  the  contact  normal,  as  a  contact 
moves  from  one  side  onto  the  next. 

DIBS  has  been  applied  to  several  analyses  of  the  effects  of  nuclear  and  chemical 
explosions  in  jointed  rocks  (Heuze,  F.  E.,  et  al.,  1990). 

B2.  CONTACTS. 

Information  on  contact  detection  and  on  comer-side  contacts  is  available  in  (Walton, 
0.  R.,  et  al.,  1991).  The  rounded-comer  contact  capability  was  developed  for  this  study. 

TV  rounded-corner  version  builds  on  the  DIBS  construct.  From  a  user’s  point  of  view 
the  model  functions  almost  the  same  as  the  sharp-comer  version.  Internally,  however,  there  are 
some  significant  differences.  Input  geometric  data  for  each  block  contain  the  coordinates  of  an 
"outer"  polygon  with  sham  comers;  this  is  just  the  polygon  the  original  DIBS  model  would  use. 
When  the  problem  is  initialized  (i.e.,  as  geometric  data  are  stored)  a  circle  of  fixed  radius,  rc  , 
is  placed  inside  each  vertex,  tangent  to  the  two  adjacent  sides.  A  new  "inner"  polygon  is 
constructed  connecting  the  centers  of  each  of  the  comer  circles.  The  side  lengths  of  the  inner 
polygon  are  exactly  the  same  as  the  lengths  of  the  line  segments  between  tangent  points  on  the 
outer  polygon  (Fig.  B-l).  The  coordinates  of  the  vertices  of  the  inner  polygon  replace  the 
original  outer  polygon  coordinates  in  the  computer  storage  arrays,  and  this  inner  polygon  is  used 
on  all  contact  searching  algorithms.  Such  a  construct  facilitates  the  use  of  most  of  the  same 
searching  algorithms  for  contact  detection  as  were  used  in  the  original  sharp-comer  version  of 


the  model.  One  major  modification  is  that  "overlap"  is  considered  to  occur  when  the  two  inner 
polygons  are  a  distance  apart  equal  to  the  sum  of  their  respective  comer  radii.  When  a  comer 
is  contacting  a  side,  the  transformed  local  coordinates,  xt  and  yf  are  still  determined  by 
(Fig.B-2): 


xt  -  (*c-*JC0S4>  +(yc-->’fl)sin<t) 


yt  -  (yc-y«)«**  -(*c-*fl)sin<i> 

The  overlap,  a,  is  now  defined  to  be 


a  -  r.t+rryt 

where  ri  and  r;-  are  the  radii  of  the  comers  for  blocks  I  and  J,  respectively  (Fig.  B-2).  This,  in 
effect,  moves  the  real  surface  (for  determination  of  contact  forces)  a  distance  rc,  outside  the 
stored  inner  polygon.  The  other  major  difference  for  the  rounded-comer  model  occurs  when  a 
contacting  comer  goes  beyond  either  end  of  the  "side"  of  the  inner  polygon,  (i.e.,  xt<0,  or 
xt>|ABj).  When  this  occurs,  the  contact  is  redefined  to  be  a  comer-comer  contact. 

The  existing  contact  forces,  including  the  currently  partially  recoverable  tangential  shear 
strain,  are  unaffected  by  the  transition  from  one  contact  type  to  another.  Only  the  indices  in  the 
linked  list  of  contact,  indicating  which  comer  and/or  side  of  which  blocks  comprise  the  contact, 
are  affected.  This  contact  information  is  in  a  separate,  dynamic,  linked  list  for  each  block.  Each 
block’s  list  contains  the  indices  of  the  sides  and/or  comers  of  blocks  contacted  by  that  block. 

B3.  CONTACT  FORCES. 

The  tangential  friction  force  model  also  was  described  by  Waiton  et  al  (  1991).  Because 
of  the  unique  character  of  the  nonlinear  normal  contact  model,  it  is  explained  in  more  detail  here. 

The  combined  effects  of  an  initially  soft,  inelastic  joint  and  the  approximate  deformation 
of  a  stiff  elastic  l  '  k  are  modeled  in  DIBS  with  a  loading  curve  that  has  a  monotonically 
increasing  slope.  The  stiffness  asymptotically  approaches  a  value  equivalent  to  the  elastic 
response  of  intact  rock  at  high  loads.  The  total  deformation  of  both  the  block  and  the  joint  is 
assigned  to  the  interface  alone,  in  the  model.  Thus,  the  apparent  overlap  between  two  blocks  in 
the  model  includes  the  displacement  associated  with  the  actual  deformation  of  the  block.  If  d 
is  the  distance  between  the  centroids  of  two  contacting  blocks,  and  Er  is  the  modulus  of  the 
blocks,  an  equivalent  spring  stiffness,  Kr  can  be  defined  between  the  centroids,  Kr  -  Erht!d  , 
where  h  is  the  block  height  and  ( is  the  joint  thickness  (see  Fig.  B-3a).  The  model  combines  this 
"rock  stiffness"  with  a  contact  relation  for  joints  that  has  a  vertical  asymptote  at  a  maximum  joint 
closure,  ac .  A  commonly  used  form  for  such  joint  behavior  (Goodman,  R.  E.,  1976  and  Bandis, 
S.  C,  et  a!.,  1983)  requires  just  two  input  parameters:  the  maximum  joint  closure,  ac  and  the 
nonnai  joint  stiffness,  Kv  The  force  displacement  assumed  for  the  joint 
is  F  *  (i^j  <xca.j)/(ac  -  a j)  ,  where  a;  is  the  displacement  (i.e.-,  closure)  of  the  joint.  The  rock 
is  assumed  to  be  linear-elastic  with  a  force  displacement  relation  F  -  Krar  ,  wh^te  ar  is  the 
displacement  associated  with  the  rock  deformation.  Rearranging  these  expressions  we  obtain. 
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3  F  +  K,a 


and  a,  -  — 

1-C  r  Kr 


for  the  displacements  of  the  joint  and  rock,  respectively.  The  total  displacement,  a,  is  simply 
the  sum  of  these  two  terms, 


Fa  p 

a  -  a;  +  a  - _ _ _ + _ 

3  r  F  +  #f1ac  Kr 
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Equation  (B-l)  may  be  solved  for  F  in  terms  of  a, 

F  -  l[aKr-ac(K,  ♦  «,)]  ♦!{[(/(:,  ♦  K,)ac-aA:r]2t4K1K,aca} 


The  slope  of  this  force  displacement  curve,  dF Id  a  ,  approaches  Kr  as  a  approaches  infinity, 

dFm^r_t _ oK,2-^  *,(*,-*!) 

rf“  2  2(a2K,2-2aaeK,(K,-fr,) » 


Figure  B-3a  is  a  schematic  representation  of  the  rock-joint  model  and  Fig.  B-3b  shows 
the  qualitative  behavior  of  the  component  and  combined  force-displacement  curves.  Various 
assumptions  are  possible  for  unloading  and  reloading  with  this  nonlinear  rock-joint  model:  (a) 
the  unload  and  reload  slopes  are  equal  to  ,  (b)  the  unload  and  reload  slopes  are  equal  to  the 
tangent  of  the  loading  curve  at  the  point  of  initial  unloading,  and  (c)  unload  and  reload  slopes 
are  set  to  a  fixed  (input)  multiplier  of  the  tangent  slope  at  initial  unloading.  The  selection  of  the 
unloading  model  determines  the  degree  of  energy  loss  during  collisions. 

B4.  OTHER  THEORETICAL  ASPECTS. 

B4.1  Viscous  Damping  Algorithm. 

Without  any  viscous  damping,  the  hysteretic  contact  model  produced  some  numerical 
noise  in  wave  propagation  calculations.  After  tests  of  various  damping  options,  we  adopted 
linear,  velocity-proportional  damping,  the  effect  of  which  was  shown  to  be  amplitude 
independent. 


F  -  Ka-DV 

The  damping  coefficient  D  is  a  user  input,  and  it  is  suggested  to  be  approximately  two  to  five 
percent  of  critical  damping 


D  -  .04  (KM)  (B-2) 

where  K  is  the  spring  stiffness  between  blocks  and  M,  the  average  block  mass.  Figure  B-4a  is 
a  plot  of  position  versus  velocity  showing  the  numerical  noise  the  code  was  generating  when  a 
square  pulse  was  applied  to  the  end  of  a  one-dimensional  chain  of  blocks,  without  any  damping. 
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Figure  B- 


In  contrast,  Figure  B-4b  is  an  illustration  of  the  same  plot  as  in  Figure  B-4a,  but  with  the 
addition  of  the  small  amount  of  viscous  damping  indicated  in  Equation  (B-2). 

B4.2  Infinite/non-reflecting  Boundaries. 

To  absorb  the  wave  with  no  reflected  pulse,  we  added  linear  damping  on  the  boundary; 
the  dashpot  coefficient  used  was  one  quarter  of  critical  damping,  as  determined  empirically 

F  -  .5  (KM)  V  (B*3)  > 

where  K  is  the  spring  stiffness  and  M  is  the  block  mass  for  blocks  next  to  the  boundary.  Note 
that  the  functional  form  of  Equation  (B-3)  is  equivalent  to  the  pC  damping  used  traditionally  in 
infinite  boundaries  for  finite  element  grids  (Lysmer,  J.  and  Kuhlemeyer,  R.  L.,  1969).  Figure  B- 
5a  is  a  plot  of  the  wave  reflected  from  the  rigid  left  boundary  of  a  one-dimensional  chain  of 
blocks,  when  the  calculation  of  Figure  B-4b  is  continued  without  the  silent  boundary  capability. 

In  contrast.  Figure  B-5b  is  the  same  plot  but  with  the  addition  of  the  silent  boundary  capability, 
and  it  shows  essentially  no  wave  reflection.  The  wave  absorption  has  proven  to  be  very  efficient 
in  two  dimensions,  as  well. 

B4.3  Other  Features.  > 

Numerous  other  features  were  put  in  DIBS  in  the  past  years.  They  are  fully  explained 
elsewhere  (Heuze,  F.  E.,  et  al.,  1990),  and  are  summarized  as: 

•  A  new  method  for  speeding  up  gravity  settling.  A  factor  of  8  in  speed  was  ^ 

obtained  by  using  a  kinetic  energy  zeroing  algorithm. 

•  A  boundary  transformation  algorithm  to  go  from  the  gravity  settling  mode  to  the 
infinite  boundary  mode. 

•  A  new  contact  search  algorithm  which  eliminates  spurious  high  forces  where  new  * 

contacts  are  discovered. 

•  A  time-history  capability  and  restart  capability. 

•  A  variety  of  mesh  generating  algorithms,  including  a  Yoronoi  polygon  generator.  > 

B5.  TIME  INTEGRATION. 

At  each  explicit  time  step  of  the  calculation,  a  loop  through  all  contacts  is  performed 
once,  accumulating  the  total  force  and  moment  on  each  particle  for  that  configuration.  Then,  the 
equations-of-motion  are  integrated  one  time  step  and  the  time  advanced.  Each  coordinate  of  each 
particle’s  position  (i.e.,  x,  y,  and  0)  is  integrated  in  time  using  a  second  order  accurate  scheme 
valid  with  non-uniform  time  steps.  For  example,  the  x-coordinate  of  the  centroid  of  a  particle 
at  time  step  n  +  1  is  determined  explicitly  from  the  positions  at  steps  n  and  n- 1  and  the 
x-direction  force  at  the  nth  time  step,  F"  ,  by  the  expression 
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where  At'  is  the  time  interval  from  n  to  b  +  1,  At  is  the  time  interval  from  n-1  to  n, 
and  i  -  F"/m  where  m  is  the  mass  of  the  block.  If  the  time  step  is  uniform  (i.e.,  At’  -  At)  , 
and  there  are  no  velocity  dependent  forces,  this  expression  reduces  to  the  familiar,  time  centered, 
third  order  accurate  Verlet  scheme  commonly  used  in  molecular  dynamics  calculations  (Allen, 
M.  P.  and  Tildesley,  D.  J.,  1987).  Similar  expressions  are  used  for  the  y-  and  0-coordinate 
integrations,  based  on  the  total  y-direction  force  and  the  total  moment  acting  on  each  block. 

Because  this  is  an  explicit  integration  scheme,  the  time  step  must  be  kept  below  a  stability 
limit  determined  by  the  stiffest  spring  and/or  the  smallest  mass  in  the  system.  For  accuracy,  the 
time  step  must  usually  be  kept  significantly  smaller  than  the  stability  limit.  Empirical  tests  of 
accuracy  have  determined  that  on  the  order  of  50  or  100  time  steps  per  physical  oscillation  cycle 
are  needed  for  the  numerically  integrated  results  to  compare  favorably  with  analytic  solutions  for 
simple  multi-block  systems  (i.e.,  ±1%). 
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DESCRIPTION  OF  FINITE  ELEMENT  METHOD  USED 
IN  UTP  BENCHMARK  ACTIVITY 
BY  CALIFORNIA  RESEARCH  AND  TECHNOLOGY 

Yoshio  Muki,  Y.  Marvin  Ito 


Cl.  INTRODUCTION. 

Finite  element  methods  which  incorporate  an  interface  logic  are  capable  of  simulating  the 
behavior  of  jointed  rock  masses  subjected  to  quasi-static  or  dynamic  loading  environments.  At 
California  Research  and  Technology  Division  (CRT)  of  The  Titan  Corporation,  analyses  of 
jointed  rock  mass  problems  for  the  UTP  Benchmark  Activity  were  performed  using 
EXCAUBUR.  Originally  developed  from  NONSAP  (Bathe,  K.  J.,  et  al.,  1974),  a  Lagrangian 
finite  element  program  for  the  solution  of  static  and  dynamic  structural  and  continuum  problems 
with  both  geometric  and  material  nonlinearity  capabilities,  EXCAUBUR  has  undergone  extensive 
development  and  validation  over  the  course  of  many  DNA  sponsored  programs. 

The  EXCALIBUR  code  uses  an  explicit  (central-difference)  time  integration  scheme 
which  allows  efficient  solution  of  ground  shock  problems  where  there  are  large  response 
gradients  in  both  time  and  space.  An  advantage  of  the  explicit  algorithm  is  the  ability  to 
incorporate  complex  constitutive  models  which  represent  detailed  inelastic  behavior  of  geologic 
materials.  The  code  also  has  a  variety  of  special  features  to  enhance  numerical  solution 
efficiency,  including  element  underintegration,  and  controls  on  the  levels  of  physical  models 
carried  by  various  parts  of  the  finite  element  grid. 

C2.  REPRESENTATION  OF  JOINTED  ROCK  MASS. 

For  analyzing  tunnel  response  in  jointed  rock  masses,  EXCALIBUR  treats  large-scale 
material  discontinuities  with  gap  formation  and  frictional  contact  using  an  interface  element  A 
recent  enhancement  incorporates  substructuring  at  the  element  level  allowing  the  code  to  capture 
micromechanical  effects  including  anisotropy  in  the  far-field  jointed  system  without  increasing 
the  global  system  of  equations.  These  methods  of  modeling  jointed  rock  masses  are  used  in  the 
UTP  benchmark  calculations.  The  first  approach,  termed  discrete  joint  modeling,  will  be 
described  next,  followed  by  a  description  of  the  second  approach  termed  composite  jointed  rock 
modeling. 

C2.1  Discrete  Joint  Modeling. 

In  the  discrete  joint  modeling  approach,  both  the  rock  mass  and  the  rock  joint  are 
represented  by  ordinary  finite  elements,  typically  four-node  quadrilaterals.  The  mechanical 
behavior  of  the  rock  joints  is  modeled  using  a  modified  version  of  the  CRT  Interface  Element 
(Ito,  Y.  M.,  et  al.,  1981).  The  interface  element  was  developed  to  treat  the  mechanical 
interactions  of  dissimilar  or  disconnected  bodies.  The  element  typically  carries  a  high  aspect 
ratio  with  a  prescribed  stiffness  in  the  direction  normal  to  the  joint  face. 
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As  specified  for  the  Benchmark  Activity,  the  normal  stiffness  of  the  rock  joint  is 
prescribed  to  increase  with  closure.  The  element  has  a  constant  elastic  shear  stiffness,  with  shear 
bearing  capacity  along  the  joint  face  following  a  non-dilatant  Coulomb-type  friction  rule 
(dependent  on  the  normal  stress).  The  element  can  be  prescribed  to  have  no  stiffness  in  tension. 
This  combination  of  features  allows  the  element  to  model  interface  stick,  slip  and  separation. 


The  constitutive  relation  governing  the  interface  element  is  a  2-by-2  matrix,  i.e.. 
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where  x  and  y  are  shear  stress  and  strain  and  aH  and  tn  are  normal  stress  and  strain.  Under 
conditions  where  the  interface  is  closed  and  no  slip  is  occurring,  the  incremental  constitutive 
relation  is 
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where  (/'and  £'  are  the  tangent  stiffnesses.  While  under  conditions  where  slip  is  occurring,  the 
relation  is 


where  ~c  is  the  cohesion  and  p  is  the  coefficient  of  friction.  And  of  course,  during  separation 
(tn  <  0),  the  stiffness  is  the  zero. 

In  this  method,  the  thickness,  location  (hence  spacing),  and  orientation  of  the  rock  joint 
are  discretely  represented  in  the  numerical  model.  11115  is  an  important  feature  of  modeling 
jointed  rock  fields  near  buried  structures  where  the  nominal  structural  dimensions  are  of  the  same 
order  as  the  nominal  joint  spacing.  Thus,  a  small  time  step  is  usually  required,  not  only  to  retain 
numerical  stability,  but  to  accurately  track  the  nonlinear  behavior  of  the  joint  system  near  the 
tunnel. 

C2.2  Composite  Jointed  Rock  Element 

In  order  to  effectively  analyze  jointed  rock  mass  in  the  far-field,  a  composite  joint/rock 
mass  element  has  been  developed  (Muki,  Y.,  et  al.,  1992).  The  composite  jointed  rock  element 
is  a  mathematical  representation  of  the  mechanical  behavior  of  a  rock  mass  with  a  regular 
distribution  of  joints.  The  element  is  intended  to  capture  the  mechanical  behavior  of  a 
characteristic  sample  of  insitu  material  comprising  a  rock  mass  interspersed  with  joint  zones. 
Each  component  within  the  composite  joint  element  is  treated  as  a  substructure  with  a  full  and 


separate  constitutive  relationship.  The  interaction  of  the  joint  zones  and  rock  mass  is  based  on 
stress  equilibrium  and  strain  compatibility  within  the  element. 

The  composite  jointed  rock  element  assumes  a  geometric  orientation  and  distribution 
between  intact  rock  and  joint  materials,  with  the  intact  material  taking  up  most  of  the  region 
occupied  by  the  element.  For  the  sake  of  this  discussion,  the  orientation  of  the  distribution  of 
joints  is  taken  to  be  horizontal  and  vertical,  although  this  restriction  can  easily  be  changed  to 
account  for  more  arbitrary  jointed  rock  mass  configurations.  The  joint  zone  represents  the  effects 
of  all  joints  with  similar  orientation  in  the  sample.  The  key  parameter  is  the  relative  thickness 
of  (distributed)  joint  material  as  compared  to  intact  rock.  Since  there  are  two  types  of  joints, 
horizontal  and  vertical,  two  characteristic  ratios  of  intact  rock  thickness  to  joint  material  thickness 
must  be  defined:  the  ratio  of  the  width  W  of  the  intact  material  to  the  width  h>  of  the  vertical 
joint  material  is  a  ■  W/w  >  >  1 ,  and  the  ratio  of  the  height  H  of  the  intact  material  to  the  height 
h  of  the  horizontal  joint  is  defined  as  p  -  H/h»  1 . 

The  superscripts  used  in  the  development  are  as  follows: 

[  /  refers  to  intact  rock  material  quantities, 

[  ]v  refers  to  vertical  joint  material  quantities, 

[  refers  to  horizontal  joint  material  quantities. 

Also,  the  symbols  Aec  refer  to  correction  strain  and  Ac j*  refer  to  out-of-balance  stress. 
Unsi.’perscripted  quantities  refer  to  overall  element  (external)  quantities  which  define  the 
composite  representation  of  the  response  of  the  various  components  of  the  jointed  rock  mass. 
These  composite  values  are  used  during  the  time  integration.  Finally,  with  the  exception  of  out- 
of-balance  stress  terms,  which  are  carried  forward  from  the  previous  time  step,  all  quantities  are 
quoted  at  the  current  time  step. 

The  following  mathematical  development  assumes  plane  strain  behavior,  with  x  denoting 
the  horizontal  direction  and  y  the  vertical  direction.  Overall  composite  strain  increments  are 
represented  by  Ae,  {Ae}  -  {Ae^  Aey  Ay^}  and  are  uniform  over  the  entire  element. 
Correction  strain  increments  are  represented  by  Aec,  {Aec}  -  {Ae^  Ae^  Ayc}.  The  correction 
strains  represent  perturbations  about  the  overall  composite  strains  ana  affect  each  component 
differently.  These  correction  strain  increments  are  used  to  adjust  the  increments  in  the  strains 
of  both  the  intact  rock  and  the  vertical  and  horizontal  joints  to  attain  internal  stress  equilibrium. 
Note  that  the  internal  strain  correction  increments  do  not  affect  the  motions  at  the  four  corner 
nodes.  Together,  seven  strain  increment  modes  (four  internal  and  three  external)  along  with  the 
two  material  distribution  parameters,  a  and  {5,  fully  describe  the  incremental  strain  state  of  each 
of  the  components. 
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The  governing  equations  for  compatibility  of  strain  increments  are 

Ae7  -  Ae,  -  —  A  el 
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a  p 


for  intact  rock. 


Ae^  -  Aex  +  Ae^ 
AE;-Ae,-iAe; 

Ayv  -  Ay  +  Ay^-Iay^ 


for  vertical  joint  material,  and 


for  horizontal  joint  material. 
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Assuming  that  changes  in  stresses  produced  from  the  strain  changes  in  the  previous  time 
step  are  not  in  equilibrium,  an  inherited  out-of-balance  in  stress  increments  Ac/  must  be  carried 


forward  and  included  in  the  calculation  for  stress  equilibrium  at  time  t.  [For  example, 

Ao‘-(a4-AoX-4,-] 


Then,  for  stress  equilibrium,  the  following  equations  must  hold 

A<  +  Ac>Ac/ 
a  Oy  -f  A  a*  -  A  a7 
At*1"  +  Atv  -  At/ 

At**  +  At*  «*  At7 

From  the  constitutive  relationship  for  intact  rock  we  have 

”  Cn  +  C- 12  ^Ey 
Ao7  -  Cj7  Ae7  +  Cjo  Ae7 
At 7  -  C33  Ay7 
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where  the  stiffness  coefficients  C?  are  assumed-isotropic  and  related  to  bulk  modulus  K*  and  shear 
modulus  Gl  by  the  usual  relations 

C,7,  *  K1  +  —  G1 
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Similar  relations  hold  for  the  vertical  joint  material  (superscript  v)  and  for  the  horizontal  joint 
material  (superscript  h). 

Combining  Equations  (C-l)  through  (C-5)  gives  a  set  of  equations  relating  intact  rock, 
vertical  joint  and  horizontal  joint  strain  increments  to  nominal  external  and  internal  strain 
increments,  a^d  corrections  in  the  stress  increments  of  the  form: 
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When  the  composite  element  is  used  in  conjunction  with  nonlinear  materials,  a 
predictor/corrector  method  (Ito,  Y,  M.,  et  al.,  1981)  is  used  to  maintain  solution  accuracy  and 
efficiency.  Briefly,  the  incremental  stiffness  relationship  in  Equation  (C-6)  is  based  on  secant 
moduli  which  are  calculated  at  the  end  of  each  time  step  for  use  in  the  equilibrium  equation  in 
the  following  time  step.  If  the  type  of  loading  changes  (c.g.,  plastic  loading/unloading  joint 
separation  closure,  a  set  of  revised  secant  moduli  can  be  calculated  for  each  of  the  components 
undergoing  the  nominal  external  strain  increment  Note  that  the  secant  moduli  are  for  an 
assumed  strain  increment  and  may  not  be  the  same  as  for  the  strain  increments  obtained  after 
solving  Equation  (C-6).  This  difference  will  lead  to  out-of-balance  stress  {  Ao6}  which  is 
added  in  the  following  time  step. 

Given  the  material  moduli,  overall  strain  increments  and  out-of-balance,  Equation  (C-6) 
can  be  solved  for  the  internal  strain  increments.  For  elastic  conditions,  Equation  (C-6)  is 
decoupled  into  two  2-by-2  systems  of  equations  which  are  easily  solved  for  {Aec}.  Insertion 
of  {Aec}  into  the  strain  compatibility  relationships  (Equation2  C-3,  C-2,  C-3)  yields  the 
predicted  equilibrium  strains  for  each  component  of  the  element.  Stresses  and  other  state 
variables  are  then  updated  by  exercising  the  appropriate  material  model  for  each  of  the 
components  in  the  new  strain  state.  By  using  the  pred;ctor/con  ector  method,  the  number  of  calls 
to  the  material  model  subroutines  can  be  minimized. 

Finally,  it  is  also  possible  and  often  desirable  to  use  the  initial  elastic  moduli  for  each  of 
the  materials  in  the  internal  equilibrium  equations  for  all  rme,  even  though  larger  approximations 
in  out-of-balance  stresses  occur,  in  order  to  reduce  the  numerical  cost  of  the  solution.  In  this 
case,  the  out-of-balance  stress  is  used  to  carry  all  nonlinear  effects  into  the  next  time  step  in  the 
integration. 

O.  BENCHMARK  CALCULATIONS. 

C3.1  Static  Calculations  (Benchmark  Problems  1,  2  and  2-S). 

Both  discrete  and  composite  joint  modeling  methods  were  used  to  analyze  three  quasi¬ 
static  Benchmark  cases.  As  will  be  shown,  except  for  some  numerical  cell  ringing,  excellent 
agreement  was  found  in  all  cases  between  the  discrete  and  composite  joint  model  calculations. 
In  all  the  quasi-static  calculations,  no  damping  was  used  in  the  equations  of  motion. 

In  Case  1  (Senseny,  P.  E.,  1991),  a  simplified  spherically  divergent  strain  path  is  used  to 
drive  both  an  intact  rock  sample  and  a  jointed  rock  sample.  For  the  jointed  rock  sample,  both 
discrete  and  composite  joint  methods  were  employed  to  obtain  results  for  comparison.  In  Case 
2  (Senseny,  P.  E.,  1991),  a  staggered  jointed  rock  array  undergoing  uniaxial  strain  loading  is 
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analyzed.  (In  this  case,  the  possibility  of  hourglass  distortion  in  a  multi-cell  grid  was  eliminated 
by  using  full  2-by-2  integration  for  the  discrete  joint  calculation.  In  all  other  cases,  one  point 
integration  was  used.)  In  the  composite  joint  calculations,  material  ratios  (a  and  p)  of  200:1 
were  used. 

Case  2-S  (Simons,  D.  E.,  1991)  was  designed  to  exercise  the  jointed  rock  mass  models 
under  combined  shear  and  compressive  loading.  In  the  discrete  model  analysis,  the  rock  material 
was  represented  using  two  constant  strain  triangles  and  the  joint  separating  the  two  was  modeled 
using  a  four  node  quadrilateral.  Although  the  original  problem  description  contained  rollers  on 
the  left  and  lower  boundaries,  in  the  present  analysis  the  boundary  conditions  were  driven  from 
all  sides  to  minimize  nonsymmetric  effects  in  the  solution  of  the  equations  of  motion.  Of  the 
twelve  independent  displacement  degrees  of  freedom  present  in  the  system,  eight  were  controlled 
by  the  specified  displacement  boundary  condition,  the  remaining  four  were  determined  by 
equations  of  motion. 

Case  2-S  was  repeated  using  the  composite  joint  model  with  very  similar  results.  The 
composite  joint  analysis  uses  only  one  element  with  eight  displacement  degrees  of  freedom.  The 
two  internal  correction  strains  are  equivalent  to  the  four  degrees  of  freedom  solved  by  the 
equations  of  motion  in  the  discrete  joint  analysis.  These  two  strain  terms  are  solved  using  the 
governing  equilibrium  equations  and  have  vastly  different  dynamic  characteristics  than  the  four 
degrees  of  freedom  solved  by  the  central  difference  method  used  in  the  discrete  joint  calculation. 
The  out-of-balance  stresses  only  appear  in  a  portion  of  the  solution  after  a  transition  (change  of 
modulus,  eluitic/plastic  loading)  and  converge  rapidly  in  typical  solutions.  The  composite  joint 
solution  under  conditions  where  there  are  no  transitions  will  give  equilibrium  exactly.  Due  to 
the  absence  of  numerical  ringing,  the  composite  joint  model  yielded  very  clean  results  with 
adequate  time  resolution,  and  reliably  determined  plasticity  and  joint  slip  in  the  solution. 

C3.2  Dynamic  Calculations  (Cases  3  and  4). 

Cases  3  and  4  are  the  field  calculations  (Senseny,  P.  E.,  1991).  In  each  case  both  discrete 
and  composite  joint  models  were  used  to  construct  a  solution.  In  the  near  field  the  joints  were 
modeled  discretely  as  specified  and  in  the  far  field  the  composite  joint  element  was  used.  All 
elements  were  evaluated  using  2-by-2  integration.  To  control  the  numerical  ringing  during 
unload,  a  time  dependent  stiffness  proportional  damping  was  used.  From  time  t  =  0  to  25 
milliseconds,  no  damping  was  used.  After  25  milliseconds  (nominally  the  TOA  of  the  peak  at 
the  tunnel  range,  damping  was  increased  linearly  in  time  to  reach  10%  of  critical  damping  at  30 
milliseconds.  The  damping  coefficient  was  held  at  10%  critical  for  the  remainder  of  the  solution 
No  lithostatic  stresses  were  assumed  to  exist  in  these  analyses. 

The  model  of  the  steel  tunnel  liner  in  Case  4  employed  a  higher-order  elasto-plastic 
element  which  simulates  the  response  in  both  thrust  and  moment.  The  element  is  an  assemblage 
of  a  membrane  component  to  match  axial  stiffness  and  a  shear/bending  component  to  match 
moment  of  inertia.  A  layer  of  interface  elements  with  the  same  normal  stiffness  behavior  used 
to  model  the  joints  but  without  shear  capacity  was  placed  between  the  tunnel  liner  and  the 
surrounding  rock  mass.  This  allows  the  structure  to  freely  slide  along  or  separate  from  the 
surrounding  jointed  rock  mass  during  deformation. 
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UNDERGROUND  TECHNOLOGY  PROGRAM  (UTP)  BENCHMARK  ACTIVITY 
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Dl.  INTRODUCTION. 

The  distinct  element  method  is  a  recognized  discontinuum  modeling  approach  for 
simulating  the  behavior  of  jointed  media  subjected  to  quasi -static  or  dynamic  conditions.  The 
method  has  three  distinguishing  features  which  makes  it  well  suited  for  discontinuum  modeling: 

(1)  The  medium  is  simulated  as  an  assemblage  of  blocks  which  interact  through 
comer  and  edge  contacts. 

(2)  Discontinuities  are  regarded  as  boundary  interactions  between  these  blocks; 
discontinuity  behavior  is  prescribed  for  these  interactions. 

(3)  The  method  utilizes  an  explicit  timestepping  (dynamic)  algorithm  which  allows 
large  displacements  and  rotations  and  general  non-linear  constitutive  behavior  for 
both  the  matrix  and  discontinuities. 

Since  the  time  the  method  was  proposed,  several  forms  of  distinct  element  codes  have  been 
developed  to  cover  a  range  of  rock  mass  strengths  and  confining  pressures  which  are  encountered 
in  situ. 


UDEC  (Universal  Distinct  Element  Code)  is  the  original  mainframe  computer  code 
developed  for  geomechanical  analysis  in  which  the  performance  of  rock  mass  may  be  dominated 
by  discontinuities  (joints,  faults,  bedding  planes).  UDEC  was  originally  developed  by  Dr.  Peter 
Cundall  through  contracts  with  the  U.S.  Army  Waterways  Experiment  Station.  Since  1983,  Itasca 
Consulting  Group,  Inc.  have  completed  a  number  of  modifications  to  the  code  which  greatly 
expand  its  range  of  applicability. 

UDEC  is  sold  by  Itasca  Consulting  Group,  Inc.  as  a  source  code  or  as  an  executable  code 
for  use  on  personal  computers  and  some  workstations.  Over  150  copies  of  the  code  have  been 
sold  world  wide.  The  most  recent  version,  Version  1.7,  was  used  as  the  basis  for  the  benchmark 
activity.  The  code  was  modified  to  include  the  intact  rock  constitutive  relation  (i.e.,  Drucker- 
Prager  model)  specified  for  the  benchmark  problem.  With  the  inclusion  of  this  constitutive 
model,  a  new  version  number  was  assigned,  Version  1.8. 

D.2  THE  DISTINCT  ELEMENT  METHOD. 

D2.1  Numerical  Formulation. 

In  the  distinct  element  method,  a  rock  mass  is  represented  as  an  assemblage  of  discrete 
blocks.  Joints  are  v. .  ,ved  as  interfaces  between  distinct  bodies.  The  contact  forces  and 
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displacements  at  the  interfaces  of  a  stressed  assembly  of  blocks  are  found  through  a  series  of 
calculations  which  trace  the  movements  of  blocks.  Movements  result  from  the  propagation 
through  the  block  system  of  a  disturbance  applied  at  the  boundary.  This  is  a  dynamic  process 
in  which  the  speed  of  propagation  is  a  function  of  the  physical  properties  of  the  discrete  system. 

The  dynamic  behavior  is  described  numerically  by  using  a  timestepping  algorithm  in 
which  the  size  of  the  timestep  is  selected  such  that  velocities  and  accelerations  can  be  assumed 
constant  within  the  timestep.  The  distinct  element  method  is  based  on  the  concept  that  the 
timestep  is  sufficiently  small  that  during  a  single  step  disturbances  cannot  propagate  from  one 
discrete  element  in  the  model  further  than  its  immediate  neighbors.  This  solution  scheme  is 
identical  to  that  used  by  the  explicit  finite  difference  method  for  continuum  numerical  analysis. 
The  timestep  restriction  applies  to  both  contacts  and  blocks.  For  deformable  blocks  (blocks 
which  are  internally  discretized  into  finite  difference  elements),  the  zone  size  is  used  to  define 
the  timestep  limitation,  and  the  stiffness  of  the  system  includes  contributions  from  both  the  intact 
rock  modulus  and  the  stiffness  at  the  contacts. 

The  calculations  performed  in  the  distinct  element  method  alternate  between  application 
of  a  force-displacement  law  at  the  contacts  and  Newton’s  second  law  of  motion  at  the  blocks. 
The  force-displacement  law  is  used  to  find  contact  forces  from  displacements.  Newton’s  second 
law  gives  the  motion  of  the  blocks  resulting  from  the  forces  acting  on  it.  If  the  blocks  are 
deformable,  motion  is  calculated  at  the  gridpoints  of  the  triangular  finite-strain  elements  within 
the  blocks.  Then,  the  application  of  the  block  material  constitutive  relations  gives  new  stresses 
within  the  elements. 

This  numerical  formulation  satisfies  momentum  and  energy  conservation  laws  by 
satisfying  Newton’s  laws  of  motion  exactly.  Although  some  error  may  be  introduced  in  the 
computer  programs  by  the  numerical  integration  process,  this  error  may  be  made  arbitrarily  small 
by  the  use  of  suitable  timesteps  and  high  precision  coordinates. 

D2.2  Rock  Joint  Representation. 

A  rock  joint  is  represented  numerically  as  a  contact  surface  (composed  of  individual  point 
contacts)  formed  between  two  block  edges.  In  genera),  for  each  pair  of  blocks  that  touch  (or  is 
separated  by  a  small  enough  gap),  data  elements  are  created  to  represent  point  contacts.  In 
UDEC,  adjacent  blocks  can  touch  along  a  common  edge  segment  or  at  discrete  points  where  a 
comer  meets  ar«  edge  or  another  comer.  For  deformable  blocks,  point  contacts  are  created  at  all 
gridpoints  located  on  the  block  edge.  Thus,  the  number  of  contact  points  can  be  increased  as  a 
function  of  the  interna!  zoning  of  the  adjacent  blocks. 

A  specific  problem  with  contact  schemes  is  the  unrealistic  response  that  can  result  when 
block  interaction  occurs  close  to  or  at  opposing  block  comers.  Numerically,  blocks  may  become 
locked  or  hung-up.  This  is  a  result  of  the  modelling  assumption  that  block  comers  are  sharp  or 
have  infinite  strength.  In  reality,  crushing  of  the  corners  would  occur  as  a  result  of  a  stress 
concentration.  Explicit  modeling  of  this  effect  is  impractical.  However,  a  realistic  representation 
can  be  achieved  by  rounding  the  comers  so  that  blocks  can  smoothly  slide  past  one  another  when 
two  opposing  comers  interact.  Comer  rounding  is  used  in  UDEC  by  specifying  a  circular  arc 
for  each  block  comer.  The  arc  is  defined  by  the  distance  from  the  true  apex  to  the  point  of 


tangency  with  the  adjoining  edges.  By  specifying  this  distance  rather  than  a  constant  radius,  the 
truncation  of  sharp  comers  is  not  severe. 

In  UDEC,  the  point  of  contact  between  a  comer  and  an  edge  is  located  at  the  intersection 
between  th“  edge  and  the  normal  taken  from  the  center  of  the  radius  of  the  circular  arc  at  the 
comer  to  the  edge.  If  two  comers  are  in  contact,  the  point  of  contact  is  the  intersection  between 
the  line  joining  the  two  opposing  centers  of  radii  and  the  circular  arcs.  The  directions  of  normal 
and  shear  force  acting  at  a  contact  are  defined  with  respect  to  the  direction  of  the  contact  normal. 
Contacts  along  the  edge  of  a  deformable  block  are  represented  by  comers  with  very  large 
rounding  lengths. 

Comer  rounding  only  applies  to  the  contact  mechanics  calculation  in  UDEC.  All  other 
calculations  and  properties  such  as  block  and  zone  mass  are  based  on  the  entire  block.  Comer 
rounding  can  introduce  inaccuracy  in  the  solution  if  the  rounding  is  too  large.  If  the  rounding 
length  is  kept  to  approximately  one  percent  (1%)  of  the  representative  block  edge  length  in  the 
model,  good  accuracy  is  achieved. 

Contact  point  information  in  UDEC  is  updated  automatically  as  block  motion  occurs.  The 
algorithms  to  perform  this  updating  must  be  computationally  efficient,  particularly  for  dynamic 
analysis,  in  which  large  displacements  may  require  deleting  and  adding  hundreds  of  contacts 
during  the  dynamic  simulation.  UDEC  takes  advantage  of  a  network  of  "domains"  created  by 
the  two-dimensional  block  assembly.  Domains  are  the  regions  of  space  between  blocks  which 
are  defined  by  the  contact  points.  During  one  timestep,  new  contacts  can  be  formed  only 
between  comers  and  edges  within  the  same  domain,  so  local  updates  can  be  executed  efficiently 
whenever  some  prescribed  measure  of  motion  is  reached  within  the  domain.  The  main 
disadvantage  of  this  scheme  is  that  it  cannot  be  used  for  very  loose  systems  because  the  domain 
structure  becomes  ill-defined. 

D2.3  Rock  Joint  Behavior. 

Numerically,  a  joint  is  a  special  contact  type  which  is  classified  as  an  edge-to-edge 
contact  in  2D.  In  UDEC,  a  joint  is  recognized  when  a  domain  is  defined  by  two  point  contacts. 
The  joint  is  assumed  to  extend  between  the  two  contacts  and  be  divided  in  half  with  each  half- 
length  supporting  its  own  contact  stress.  Incremental  normal  and  shear  displacements  are 
calculated  for  each  point  contact  and  associated  length. 

UDEC  has  two  joint  behavior  relations  to  describe  the  mechanical  response  at  the 
interface.  The  basic  joint  model  (Coulomb  slip  model)  captures  several  of  the  features  which  are 
representative  of  the  physical  response  of  joints.  In  the  normal  direction,  the  stress-displacement 
relation  is  assumed  to  be  linear  and  governed  by  the  stiffness  kn  such  that 


where  an  is  the  effective  normal  stress,  and  un  is  the  normal  displacement. 

There  is  also  a  limiting  tensile  strength,  T,  for  the  joint.  If  the  tensile  strength  is  exceeded 
(i.e.,  if  on  <  -7),  then  an  =  0.  Similarly,  in  shear  the  response  is  controlled  by  a  constant  shear 


stiffness,  ks .  The  shear  stress,  xs ,  is  limited  by  a  combination  of  cohesive  (C)  and  frictional  (4>) 
strength.  This,  if 

\Tt\sC*onm+~Tam 
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then 
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or  else,  if 
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where  u*  is  the 
displacement. 

elastic  component  of  the  shear  displacement,  and  us 

is  the  total  shear 

For  the  benchmark  activity,  a  more  comprehensive  joint  model,  called  the  continuously- 
yielding  joint  model  (Cundall,  P.  A.  and  Lemos,  J.  V.,  1990)  was  used.  The  continuously- 
yielding  joint  model,  proposed  by  Cundall  and  Hart  (1984)  and  revised  by  Lemos  (1987),  is 
intended  to  simulate  the  intrinsic  mechanism  of  progressive  damage  of  the  joint  under  shear. 
This  approach  produces  consistent  responses  in  the  varied  conditions  encountered  in  numerical 
modeling.  The  model  also  provides  continuous  hysteretic  damping  for  dynamic  simulations. 

The  response  to  normal  loading  is  expressed  incrementally  as 

A  cl  -  kn  A u„ 

FI  FI  FI 

€n 

where  kn  is  the  normal  stiffness,  given  by  kn  -  a.non  ,  a  simple  relation  representing  the 
observed  increase  of  stiffness  with  normal  stress,  where  an  and  en  are  model  parameters.  The 
user  may  specify  maximum  and  minimum  normal  stiffnesses  and,  in  general,  zero  tensile  strength 
is  assumed. 

One  feature  of  the  continuously-yielding  model  is  the  ability  to  simulate  the  intrinsic 
mechanisms  of  progressive  damage  of  the  joint  under  shear.  Since  this  feature  was  not  required 
for  the  benchmark  activity,  the  shear  stress-displacement  relation  previously  described 
(Equations  D-l  through  D-4)  was  used. 

D2.4  Block  Deformability. 

In  UDEC,  each  block  is  automatically  discretized  into  triangular  constant-strain  elements. 
These  elements  may  follow  an  arbitrary,  nonlinear  constitutive  law  (e.g.,  Mohr-Coulomb  failure 
criterion  with  non-associated  flow  rule).  Other  nonlinear  plasticity  models  recently  added  to 
UDEC  include  a  ubiquitous  joint  model  and  strain-softening  models  for  both  shear  and 
volumetric  (collapse)  yield.  As  mentioned  in  the  introduction,  the  Drucker-Prager  model  was 


2ddcd  specifically  for  use  in  the  benchmark  activity.  The  complexity  of  deformation  of  the 
blocks  depends  on  the  number  of  elements  into  which  the  blocks  are  divided.  For  the  benchmark 
activity,  either  4  or  16  triangular  zones  were  used  for  each  square  block. 

D2.5  Numerical  Damping. 

If  natural  energy  dissipation  such  as  inter-block  sliding  or  internal  block  failure 
accompanies  the  discontinuum  analysis,  unwanted  vibrations  due  to  initial  or  transient  force 
imbalance  will  be  absorbed.  However,  if  the  analysis  is  predominantly  elastic,  some  artificial 
damping  will  be  necessary.  Damping  is  used  in  UDEC  to  solve  both  static  and  dynamic 
problems.  Static  problems  generally  require  more  damping  than  dynamic  ones.  Two  types  of 
damping,  mass-proportional  and  stiffness-proportional,  are  available  in  UDEC.  Mass-proportional 
damping  (or  viscous  damping)  applies  to  centroids  of  rigid  blocks  or  gridpoints  of  deformable 
blocks  a  force  which  is  proportional  to  the  (mass)  velocity  but  in  the  opposite  direction. 
Stiffness-proportional  damping  applies  to  contacts  or  stresses  in  zones  as  a  force  which  is 
proportional  to  the  incremental  force  or  stress  and  in  the  same  sense.  Either  form  of  damping 
may  be  used  separately  or  in  combination.  The  use  of  both  forms  of  damping  in  combination 
is  termed  Rayleigh  damping  (Bathe,  KL  J.  and  Wilson,  E.  L.,  1976). 

Qualitatively,  the  mass-proportional  part  tends  to  act  on  lower  frequency  modes  which 
usually  are  associated  with  the  movement  in  unison  of  several  blocks  or  gridpoints  ("sloshing”), 
while  the  stiffness-proportional  component  damps  higher  frequency  inter-block  vibrations 
("rattling"). 

The  mass  damping  force  terms  {d},  for  translational  degrees  of  freedom  in  the  momentum 
equation  take  the  form 

{d}  - 

where  \M]  is  the  mass  matrix,  {«}  is  the  velocity  vector,  and  a  is  the  mass-proportional 
damping  constant. 

Vibrational  energy  generated  at  contacts  between  blocks  or  within  deformable  zones  is 
damped  by  applying  stiffness-proportional  damping  at  the  contact  or  zone.  The  damping  force 
is 

{s}  -  p  [*]{*} 

where  J7C]  is  the  stiffness  matrix,  and  p  is  the  stiffness-proportional  damping  constant.  The 
stiffness  damping  force  is  omitted  if  sliding  occurs  at  the  contact  or  if  failure  occurs  within  the 
zone  because  frictional  dissipation  provides  natural  damping. 

For  a  multiple  degree-of-ffeedom  system,  choice  of  the  damping  constants,  a  and  P, 
cannot  be  made  with  certainty.  However,  the  critical  damping  ratio,  y,  at  any  natural  (angular) 
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frequency  of  the  system,  o>,  can  be  found  (Bathe,  K.  J.  and  Wilson,  E.  L.,  1976)  to  be 

Y  -  A  —  +  pto 
2[o) 

The  level  of  damping  is  seen  to  be  frequency-dependent  The  values  of  a  and  p  must  be  chosen 
to  provide  a  suitable  fraction  of  critical  damping. 


Equation  (D-5)  reaches  a  minimum  at 
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The  fundamental  frequency  is  then  defined  as 

■/min  “ 


where  units  are  cycles/second.  The  values  of  Ymjn  and  /mjn  are  required  input  for  damping  in 
UDEC.  For  static  analysis,  the  values  have  often  been  determined  either  by  using  a  simplified 
analog  of  the  model  (for  example,  the  vibration  of  an  equivalent  elastic  half  space)  or  by 
monitoring  a  short  undamped  run  so  that  dominant  modes  for  damping  can  be  identified.  An 
alternative  approach  used  in  UDEC  is  to  employ  an  "adaptive"  damping  scheme  which  adjusts 
the  mass  damping  constant,  a,  automatically  to  the  changing  conditions  of  the  problem  during 
solution  (Cundall,  P.  A.,  1982).  The  adaptive  damping  algorithm  monitors  the  rate  of  energy 
change  in  the  model  during  solution.  The  value  of  a  is  adjusted  as  a  function  of  the  rate  of 
energy  dissipation  and  the  rate  of  change  of  kinetic  energy  in  the  model.  The  algorithm  is 
similar  to  a  servo-mechanism  in  which  an  output  system  control  parameter  is  continuously 
adjusted  based  on  the  measured  system  response. 


For  dynamic  analysis,  the  damping  should  attempt  to  reproduce  the  frequency-dependent 
damping  of  natural  materials  at  the  correct  levels.  For  geologic  materials,  this  is  generally  2% 
to  5%  of  the  critical  damping.  The  dominant  frequencies  are  a  combination  of  the  input  wave 
frequencies  and  the  natural  modes  of  the  system.  For  the  benchmark  activity,  no  numerical 
damping  was  used  except  to  arrive  at  the  prescribed  static  in-situ  conditions  prior  to  application 
of  the  dynamic  load. 


D2.6  Solution  Stability. 

As  mentioned  previously,  the  solution  scheme  used  for  the  distinct  element  method  is  only 
conditionally  stable.  The  limiting  timestep  is  determined  which  satisfies  both  the  stability 
criterion  for  calculation  of  internal  block  deformation  as  well  as  that  for  the  inter-block  relative 
displacement.  The  approximate  timestep  required  for  the  stability  of  block  deformation 
computations  is 

A  tn  -  2  min 
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where  mt  is  the  mass  associated  with  block  node  i,  and  kt  is  the  measure  of  stiffness  of  the 
elements  surrounding  the  node.  The  ratio  of  mass  to  stiffness  is  analogous  to  the  highest 
eigenfrequency,  <omax ,  of  a  linear  elastic  system. 


The  stiffness  term  kt  must  account  for  both  the  stiffness  of  the  intact  rock  and  that  of  the 
discontinuities.  It  is  calculated  as  the  sum  of  the  two  components: 


The  first  term  on  the  right-hand  side  represents  the  sum  of  the  contributions  of  the  stiffness  of 
all  elements  connected  to  node  i,  which  are  estimated  as 


^  3 


k+1g 


1  b. 


max 


min 


where  K  and  G  are  the  bulk  and  shear  elastic  moduli  of  the  block  material,  respectively,  6max  is 
the  largest  zone  edge,  and  /tmin  is  the  minimum  height  of  the  triangular  element.  The  joint 
stiffness  term  kyt  exists  only  for  nodes  located  on  the  block  boundary,  and  is  taken  as  the 
product  of  the  normal  or  shear  joint  stiffness  (whichever  is  larger)  and  the  sum  of  the  lengths  of 
the  two  block  edge  segments  adjacent  to  node  i. 


For  calculations  of  inter-block  relative  displacement,  the  limiting  timestep  is  calculated, 
by  analogy  to  a  simple  degree -of-freedom  system,  as 


Atb  «  (frac)  2 


where  Afmjn  is  the  mass  of  the  smallest  block  in  the  system,  and  12^  is  the  maximum  contact 
stiffness.  The  factor  frac  is  a  user-supplied  value  which  accounts  for  the  fact  that  a  single  block 
is  in  contact  with  several  blocks.  A  typical  value  for  frac  is  0.1. 

The  controlling  timestep  for  a  distinct  element  analysis  is 

At  =  min(Atn,Atb) 

When  stiffness-proportional  damping  is  used,  this  timestep  restriction  may  not  guarantee  stability. 
Belytschko  (1983)  shows  that  the  timestep  in  this  case  should  be  adjusted  as  follows: 

At  =  — - — [l/l  +Y2  ~y] 


where  o)max  is  the  highest  eigenfrequency  of  the  system,  and  y  is  the  fraction  of  damping  at  this 
frequency.  This  is  accounted  for  in  the  UDEC  analysis  by  a  user-supplied  factor  based  on  the 
selected  damping  conditions. 
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APPENDIX  E 


THE  APPLICATION  OF  THE  FLEX  PROGRAM 
TO  UTP  BENCHMARK  PROBLEMS 

Felix  Wong,  David  Vaughan 


El.  INTRODUCTION. 

Weidiinger  Associates  analyzed  the  series  of  UTP  Benchmark  Problems  using  the  FLEX 
finite  element  program  (Vaughan,  D.  K.  and  Richardson,  E.,  1991).  The  FLEX  code  was 
developed  by  Weidiinger  Associates  to  solve  a  wide  range  of  structural  response,  structure-media 
interaction  and  continuum  wave  propagation  problems. 

Several  of  the  smaller  benchmark  problems  were  solved  running  FLEX  on  an  IBM 
RS/6000  workstation.  The  larger  models  were  analyzed  on  a  supercomputer  at  Los  Alamos 
National  Laboratory. 

E2.  OVERVIEW  OF  THE  FLEX  CODE. 

FLEX  is  a  finite  element  code  designed  for  2-D  and  3-D  dynamic  analysis  using  an 
explicit  time  integration  approach.  Both  geometric  and  material  nonlinearity  may  be  included 
in  an  analysis.  Static  problems  may  also  be  solved  using  a  dynamic  relaxation  approach. 

E2.1  Time  Integration  Approach. 

The  equilibrium  equations  of  motion  for  a  finite  element  grid  can  be  expressed  as: 

[M]u  +  [C]u  +  [K]u  -  Q  (E-1) 

where  u  is  the  displacement  vector  for  all  nodes  of  the  grid  (the  single  and  double  dots  indicate 
single  and  double  time  derivatives,  respectively.  [Af|,  [C]  and  [X]  are  the  mass,  damping  and 
stiffness  matrices  and  Q  is  a  vector  of  applied  nodal  forces.  FLEX  solves  Equation  (E-l)  using 
the  second  order,  central  difference  approach  to  explicitly  integrate  the  equations  of  motion 
forward  in  time.  A  diagonal  lumped  mass  matrix  is  used  and  a  diagonalized  stiffness  and/or 
mass  proportional  damping  matrix  is  also  assumed.  This  approach  avoids  the  necessity  of 
assembling  and  manipulating  large  global  stiffness,  mass  and  damping  matrices  and  provides  a 
computationally  efficient  approach  to  the  solution  of  large,  nonlinear  wave  propagation  problems. 

The  adopted  approach  is  numerically  stable  for  time  steps  which  are  no  greater  than  the 
Courant  stability  limit.  This  criteria  can  be  thought  of  as  the  time  required  for  a  wave  to 
propagate  across  the  minimum  dimension  of  an  element  within  the  model.  Since  a  model  may 
contain  locally  small  or  stiff  elements  which  tend  to  reduce  the  allowable  time  step  for  an 
element,  an  approach,  termed  subcycling,  is  used  to  allow  different  time  steps  in  different  regions 
of  the  model.  This  approach  minimizes  the  time  step  penalty  associated  with  locally  small  or 
stiff  elements. 
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E2.2  Element  Library. 


FLEX  may  be  used  to  solve  both  2-D  and  3-D  models.  For  2-D  plain  strain  or 
axisymmetric  models,  the  available  elements  types  are:  4  node  quadrilateral  continuum  elements, 
2  node  layered  shell  (Mindlin  plate)  and  2  node  spring  and  dashpot  continuum  elements.  For  3-D 
models,  the  available  elements  types  are:  8  node  hexahedron  continuum  elements,  4  node  layered 
shell  and  2  node  spring  and  dashpot  elements. 

All  elements  assume  linear  shape  functions  and  utilize  single  point  quadrature  integration 
for  all  integrals  over  the  volume  of  an  element.  Since  this  approach  under  integrates  the  internal 
strain  energy  of  an  element,  there  are  deformational  modes,  termed  hourglass  modes,  of  an 
element  which  are  unresisted.  The  adopted  formulation  for  the  continuum  elements  (Flanagan, 
D.  P.  and  Belytschko,  T.,  1981)  includes  orthogonal  hourglass  control  to  suppress  these  hourglass 
modes. 

E2.3  Constitutive  Models. 

A  number  of  constitutive  models  are  available  within  FLEX  to  allow  the  simulation  of 
a  material’s  nonlinear  constitutive  behavior  in  a  high  stress  environment.  These  models, 
developed  by  Weidlinger  Associates,  include  the  cap  family  of  models  for  soil  and  rock 
(Dimagio,  F.  L  and  Sandler,  I.  S.,  1971)  and  equivalent  models  for  concrete  (Levine,  H.  S., 
1982).  These  constitutive  models  use  an  incremental  plasticity  formulation  with  a  fixed  failure 
surface  and  moveable  cap.  An  associated  flow  rule  is  assumed.  This  approach  provides 
guarantees  of  uniqueness,  continuity  and  stability  for  material  behavior  when  subjected  to 
multi-dimensional  strain  paths.  Both  rate  dependent  and  rate  independent  models  are  available. 
Also,  two-  and  three-invariant  plasticity  surfaces  are  available. 

E2.4  Geometric  Nonlinearity. 

Beth  small  rotation/small  strain  and  large  rotation/moderate  strain  options  are  available 
within  the  program.  The  large  rotation/moderate  strain  option  uses  an  updated  Lagrangian 
approach  to  account  for  the  geometric  nonlinear  behavior  of  the  model. 

E2.5  Joint  Modeling. 

FLEX  provides  two  methods  (Levine,  H.  S.,  et  al.,  1989)  for  modeling  joints.  These  are 
the  continuum  interface  element  approach  and  the  slide  line  approach. 

The  interface  element  approach  uses  ’thin’  continuum  elements  to  represent  a  joint 
between  two  blocks  of  material  in  an  analysis.  The  interface  elements  are  not  the  actual 
thickness  of  the  joint  since  this  would  ;esult  in  severe  time  step  penalties  and  aspect  ratio 
problems.  Instead,  an  interface  element  simulates  a  region  of  continuum  material  with  a  joint 
embedded  within  it.  These  elements  have  special  joint  constitutive  relations  in  addition  to  the 
standard  constitutive  properties  of  the  continuum  material  on  either  side  of  the  joint.  Simulation 
of  slip  and  separation/recontact  can  be  achieved  using  this  approach.  The  normal  and  shear 
stresses  of  interface  elements  are  modified  based  upon  whether  contact  or  separation  is  detected 
and  by  limiting  the  shear  stresses  tangential  to  the  boundary  using  a  Mohr-Coulomb  failure 
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surface  that  is  fit  to  joint  data.  Because  the  nodal  connectivity  of  interface  elements  do  not 
change  during  an  analysis,  this  approach  is  inherently  restricted  to  problems  in  which  the  amount 
of  sliding  along  an  interface  is  much  less  than  one  element  dimension. 

The  slide  line  approach  used  in  FLEX  is  a  penalty  function  approach  (Halquist,  J.  O.,  et 
al.,  1985).  This  approach  allows  the  user  to  define  multiple  domains  within  the  computational 
model  which  interact  with  each  other  through  the  forces  produced  by  contact  (penalty)  springs. 
Penalty  springs  form  between  two  domains  whenever  the  nodes  of  one  domain  penetrate  the 
exterior  boundary  of  another  domain.  This  approach  can  be  used  for  problems  which  include 
significant  sliding  along  a  joint  since  it  accounts  for  an  adaptive  connectivity  which  changes  as 
nodes  slide  beyond  the  current  contact  element  on  the  boundary  of  a  domain. 

The  main  features  of  the  slide  line  contact  computation  are  as  follows:  At  each  time  step, 
after  the  geometry  of  the  slide  line  boundary  nodes  have  been  updated  from  the  last  time  step: 

(1)  The  closest  element  on  the  adjacent  domain  is  determined  for  each  slide  line  node. 

(2)  For  slide  line  nodes  which  are  in  contact  with  the  adjacent  domain:  the  normal 
penetration,  e,  into  the  adjacent  domain  and  the  increment  of  sliding,  6,  along  the 
adjacent  domain’s  perimeter  during  the  present  time  step  is  computed. 

(3)  The  normal  contact  force,  F  ,  due  to  a  nodes  penetration  of  an  adjacent  domain 
at  time  n  is  computed  as 

F  n  -  kznA 

where  k  represents  the  normal  stiffness  of  the  joint  and  A  is  the  tributary  contact 
area  associated  with  the  node.  Currently,  k  is  considered  to  be  a  constant  but  this 
could  be  modified  so  that  k  =  /(e). 

(4)  Equal  and  opposite  reaction  forces  are  distributed  to  the  nodes  which  define  the 
penetrated  element  of  the  adjacent  domain.  Once  this  process  is  complete  for  all 
nodes  on  both  sides  of  the  slide  line,  the  total  normal  interaction  force,  /,  at  a 
node  is  known.  /  includes  F  for  the  node  and  any  reaction  forces  due  to  nodal 
penetration  from  the  adjacent  domain. 

(5)  The  normal  contact  stress,  o,  for  each  node  at  time  n  is  defined  as: 


A 


(6)  The  joint  shear  stress,  x,  at  time  n  is  evaluated  incrementally  as: 

x*  ■>  x”-1  +  g 6 
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where  g  is  the  shear  stiffness  of  the  joint,  assumed  to  be  a  constant,  x  must  also 
satisfy  the  Coulomb  friction  limit  defined  as: 

|x”|  sC  +  c^taml* 

where  C  is  the  initial  cohesion  and  /  is  the  shear  friction  angle  for  the  joint. 

(7)  Nodal  forces  due  to  x  are  computed  for  each  node  and  equal  and  opposite  reaction 
forces  are  distributed  to  the  nodes  of  the  adjacent  contact  surface. 

(8)  Slide  line  processing  is  complete  for  this  time  step. 

Although  the  amount  of  joint  slip  anticipated  to  occur  for  the  benchmark  problems  was 
not  great,  the  slide  line  approach  was  chosen  to  be  used  in  the  UTP  Benchmark  series  since  it 
was  anticipated  that  future  problems  of  interest  to  DNA  may  include  large  amounts  of  joint  slip 
which  would  be  outside  the  range  of  applicability  of  interface  element  modeling.  Consequently, 
the  benchmark  problems  were  seen  as  an  opportunity  to  gain  practical  experience  modeling  joints 
using  slide  lines. 

E2.6  Artificial  Viscosity. 

FLEX  allows  the  use  of  either  linear  or  quadratic  artificial  viscosity  for  shock  related 
wave  propagation.  These  viscosities  allow  the  smoothing  out  of  shocks  which  form  in  an 
analysis  due  to  the  hardening  behavior  of  certain  types  of  materials  such  as  many  clays  .  These 
controls  provide  a  way  to  control  spurious  numerical  oscillations  within  the  model  which  occur 
due  to  the  fundamental  frequency  limitation  of  computational  models.  As  is  typical  for  most 
codes,  the  viscosity  is  applied  to  the  volumetric  strain  rate  and  is  used  for  compressive  strain 
rates  only. 

E2.7  Boundary  Conditions. 

A  number  of  boundary  condition  options  are  available  within  the  program  including 
prescribed  pressure,  prescribed  velocity  and  absorbing  boundary  conditions.  The  absorbing 
boundary  condition  uses  a  standard  normally  incident  plane  wave  assumption  (Lysmer,  J.  and 
Kuhlemeyer,  R.  L.,  1969).  The  methodology  is  derived  from  one  dimensional  wave  propagation 
impulse-momentum  balance  relations.  The  applicable  relationship  for  the  normal  velocity  at  an 
absorbing  boundary  is: 
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where  Vn  is  the  normal  velocity,  a  is  the  normal  stress,  p  is  the  mass  density  of  the  material  and 
cp  is  the  dilatational  wave  speed  of  the  material.  The  equivalent  relationship  for  the  tangential 
motion  at  the  boundary  is: 
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where  Vt  is  the  tangential  velocity,  x  is  the  tangential  shear  stress,  and  cs  is  the  shear  wave  speed 
of  the  material. 

E3.  APPLICATION  OF  FLEX  TO  BENCHMARK  PROBLEMS. 

This  section  describes  specific  modeling  issues  and  approaches  which  were  used  for  the 
series  of  benchmark  problems.  As  mentioned  previously,  the  slide  line  approach  was  used  to 
explicitly  model  joint  behavior  in  all  the  benchmark  analyses.  No  stiffness  or  mass  proportional 
damping  was  used.  Standard  values  of  artificial  viscosity  were  used.  All  analysis  were 
performed  using  2-D  plane  strain  or  axisymmetric  models. 

The  first  three  benchmark  problems  are  all  composed  of  only  a  few  blocks  with  explicitly 
defined  joints  separating  them.  The  last  two  benchmarks  involved  simulations  which  included 
regions  of  the  model  with  explicitly  modeled  blocks  and  joints  and  other  regions  with  implicitly 
jointed  material. 

Modeling  problem  2  produced  excessive  numerical  chatter  in  the  joints.  Improvements 
in  understanding  the  interaction  of  discrete  blocks  using  slide  line  joint  models  allowed  for 
improved  modeling  of  the  final  two  benchmark  problems  resulting  in  negligible  numerical  noise 
and  chatter  due  to  joint  contact. 

E3.1  Constitutive  Model. 

The  requested  Drucker  Prager  constitutive  model  to  be  used  for  the  benchmark  problems 
was  represented  within  the  program  by  fitting  the  failure  surface  of  the  standard  two  invariant 
cap  model  to  the  Drucker  Prager  failure  surface.  The  cap  portion  of  the  model  was  deactivated 
for  all  calculations.  The  standard  associated  flow  rule  was  utilized. 


E3.2  Explicitly  Jointed  Regions. 

In  regions  of  the  model  where  blocks  and  their  surrounding  joints  were  explicitly 
modeled,  all  blocks  were  modeled  as  fully  deformable  solids  with  an  8  x  8  arrangement  of  4 
node  quadrilateral  continuum  elements.  Slide  lines  were  placed  on  all  sides  of  each  block  to 
simulate  the  joint  behavior.  The  normal  stiffness,  k,  and  shearing  stiffness,  g,  of  the  joint  were 
assumed  constant  Although  the  slide  line  penalty  springs  could  have  been  altered  to  model  the 
specified  nonlinear  stiffness  of  the  joints,  the  linear  approximation  was  considered  adequate  and 
facilitated  matching  the  overall  stiffness  of  the  implicit  and  explicitly  jointed  regions  throughout 
the  range  of  applied  pressure  loading. 

E3.3  Implicitly  Jointed  Regions. 

The  implicitly  jointed  material  accounts  for  the  presence  of  joints  within  a  continuum 
region  of  the  model  by  altering  the  constitutive  parameters  of  the  native  continuum  material. 
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This  approach  (Si%;b.  &,,  1973)  alters  the  elastic  stiffness  properties  of  the  continuum  so  that 
they  represent  the  composite  behavior  of  multiple  layers  of  intact  material  separated  by  joints 
with  independent  stiffness  characteristics.  The  modified  continuum  properties  are  assumed  to 
remain  isotropic,  i.e.  the  continuum  region  is  assumed  uniformly  jointed  in  all  directions.  The 
elastic  bulk  and  shear  modulus  of  the  intact  continuum  material  was  K=  20  GPa  and  G=  12  GPa. 
For  implicitly  jointed  material,  me  reduced  elastic  properties  were  K1  =  9.5  GPa 
and  G/  =  8.36  GPa.  The  equivalence  of  the  implicit  and  explicitly  jointed  material  was 
demonstrated  by  performing  1-D  rock  column  simulations  for  both  representations  and  comparing 
time  histories  of  stress  and  velocities  at  equivalent  spatial  locations. 

E3.4  Static  Overburden  Solution. 

A  static  solution  of  the  model  subjected  to  gravitational  loading  was  performed  to  provide 
the  insitu  stress  state  as  initial  conditions  for  the  dynamic  response  analysis  of  Benchmark 
Problems  3  and  4.  The  geometrical  spreading  included  within  the  dynamic  model  was  a 
complication  in  performing  a  proper  static  analysis.  This  difficulty  was  avoided  by  wrapping 
additional  layers  of  elements  around  the  dynamic  model  in  order  to  p  roduce  a  static  model  with 
a  globally  oriented,  rectangular  outer  perimeter.  Tins  model,  subjected  to  gravitational  load,  was 
then  solved  using  dynamic  relaxation.  The  dynamic  analysis  used  the  stress  distribution  from 
the  static  analysis  as  appropriate  initial  conditions. 

The  dynamic  response  of  the  two  benchmark  problems  with  and  without  insitu  initial 
conditions  was  investigated.  It  was  found  that  the  insitu  stress  state  was  not  very  important  for 
this  group  of  problems  due  to  the  relatively  low  stress  levels  compared  ?o  the  dynamic  loading 
applied. 

E3.5  Boundary  Conditions. 

Standard  boundary  conditions  (prescribed  pressure  time  histories,  symmetry  p'anes,  etc.) 
were  used  for  all  benchmark  problems.  One  item  of  note  is  the  proper  use  of  the  absorbing 
boundary  at  the  bottom  of  the  model  for  benchmark  problems  4  and  5.  The  default  diiaiational 
wavespeed  used  for  the  absorbing  conditions  (Equations  E-2  and  E-3)  are  the  elastic  wave  speeds 
of  the  materia!.  Since  the  loading  wave  speed  of  the  continuum  material  for  the  applied  loading 
is  significantly  slower  than  the  elastic  wave  speed,  the  bottom  boundaries  tended  to  reflect  more 
of  the  applied  pressure  wave  than  was  desirable.  Consequently,  a  proper  loading  wave  speed  was 
determined  for  the  material  to  more  appropriately  match  the  impedance  characteristics  of  the 
absorber  to  that  of  the  discretized  domain.  Standard  program  options  allow  the  user  to  override 
the  default  elastic  wavespeeds  for  absorbing  boundary  conditions  when  the  need  arises. 
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E4.  LESSONS  LEARNED  UNDER  UTP  BENCHMARK  SERIES  FOR  MODELING 

JOINTED  BLOCK  BEHAVIOR.' 

The  UTP  series  of  benchmark  problems  provided  a  useful  set  of  simulation  problems  and  1 

a  useful  environment  for  evaluating  and  enhancing  joint  modeling  simulation  capabilities.  The 
benchmark  series  consisted  of  a  group  of  problems  ranging  from  small,  single  element  problems 
to  large  jointed  block  assemblages.  The  solutions  for  the  simpler  problems  could  be  evaluated 
analytically  to  validate  the  numerical  solutions.  The  more  complex  problems  could  be  checked 
for  rational  consistency  and  compared  with  other  independent  numerical  solutions.  * 

While  progressing  through  the  series  of  benchmark  problems,  WA  gained  some  useful 
insight  into  the  application  of  its  simulation  software  to  jointed  rock  problems.  This  write  up 
briefly  describes  the  primary  lessons  learned  during  the  benchmark  study  phase. 

» 

E4.1  Representation  of  Joint  Stiffness. 

The  penalty  function  slide  line  approach  used  in  FLEX  uses  a  default  penalty  spring 
stiffness  related  to  the  geometric  size  of  the  elements  adjacent  to  the  slide  line  and  their  elastic 
stiffness  properties.  This  provides  an  internally  computed  spring  stiffness  which  minimizes  the 
amount  of  slide  line  penetration  which  may  occur  but  still  provides  a  reasonable,  stable  time  step 
for  the  calculation.  Initial  benchmark  calculations  using  default  penalty  spring  values  showed 
them  to  be  unacceptable  when  trying  to  simulate  the  macro  effective  stiffness  of  a  jointed  block 
region  when  the  prescribed  joint  stiffness  differs  significantly  from  the  default  value. 

The  code  input  was  generalized  to  allow  a  user  prescribed  linear  penalty  spring  stiffness,  > 

K.  This  allows  a  user  to  match  available  joint  stiffness  data.  This  provides  for  a  much  better 
representation  of  the  stiffness  characteristics  of  the  jointed  block  assemblage. 

E4  2  Numerical  Chatter  at  Joints. 

> 

Our  initial  simulation  of  benchmark  problem  #2  manifested  a  significant  amount  of 
numerical  chatter  or  vibration.  Numerical  experiments  using  damping  to  control  this  chatter  were 
unsatisfactory.  The  problem  was  closely  scrutinized  to  identify  and  understand  the  source  of  the 
numerical  noise.  It  was  determined  that  there  were  two  contributing  issues: 

(1)  Numerical  precision 

(2)  Treatment  of  slide  line  joints  at  block  comers 
E4.2.1  Numerical  Precision  Issues  Related  to  Joint  Modeling. 

S 

A  recently  acquired  IBM  RS/6000  workstation  was  used  to  perform  the  small  benchmark 
problem  simulations.  The  larger  benchmark  problems  were  simulated  on  a  CRAY  computer  at 


lrrhis  subsection  was  provided  during  the  review  of  the  initial  draft  of  the  report  to  clarify 
certain  modifications  implemented  during  the  course  of  the  study. 
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> 


Los  Alamos.  Benchmark  #2  was  run  on  the  IBM  using  a  single  precision  version  of  the  FLEX 
code. 


We  determined  that  for  certain  classes  of  problems,  using  a  32  bit  precision  algorithm  for 
slide  line  computations  results  in  a  significantly  noisier  response  than  if  64  bit  arithmetic  were 
used.  This  results  from  the  fact  that  the  amount  of  penetration,  e,  which  is  computed  by  the  slide 
line  logic  is  based  on  the  global  coordinates  of  slide  line  nodes.  Small  amounts  of  penetration 
which  would  result  in  the  initiation  of  resisting  penalty  function  forces  in  a  64  bit  computation 
can  result  in  an  evaluation  of  no  penetration  using  32  bit  logic.  Consequently,  64  bit 
computations  tend  to  have  a  smoother  contact  initiation  phase  and  a  correspondingly  smoother 
response. 

Follow-on  computations  were  performed  using  64  bit  precision  and  resulted  in  reduced 
noise  and  improved  slide  line  response. 

E4.2.2  Joint  Comer  Treatment  for  Blocks. 

Another  contribution  to  the  numerical  chatter  observed  in  Bench  ^rk  #2  was  related  to 
the  treatment  of  block  comers.  The  jointed  rock  benchmark  problems  contain  a  number  of 
discrete  blocks,  each  with  joint  planes  all  around.  Each  block  is  composed  of  a  very  stiff 
material,  which  has  the  potential  to  exhibit  a  fair  amount  of  elastic  ringing.  The  basic  penalty 
function  slide  line  approach  defines  the  penalty  spring  forces  acting  on  a  node  based  on  the 
amount  of  the  node’s  penetration  into  the  surface  of  an  adjacent  domain  or  block.  When  a  node 
slides  beyond  the  extent  of  an  adjacent  block’s  comer,  it  no  longer  can  penetrate  that  block  and 
therefore  it  is  no  longer  in  contact  with  the  block.  Consequently,  any  previous  penalty  spring 
forces  which  were  applied  to  the  node  are  removed.  This  sudden  transition  from  a  penetrating 
node  to  a  non-penetrating  node  can  produce  a  local  ringing  of  the  node.  Since  the  node  is  in 
close  proximity  to  the  adjacent  block’s  comer,  the  local  ringing  of  the  node  can  move  it  back  into 
contact  with  the  block,  causing  an  amplification  of  the  numerical  oscillation  of  the  node.  The 
larger  the  number  of  block  comers  within  the  model,  the  larger  the  number  of  potential  sources 
of  numerical  noise  from  this  effect. 

Upon  careful  considerations  of  the  fundamental  physical  phenomenon  we  are  trying  to 
simulate  for  the  jointed  block  problem,  it  was  obsep/ed  that  the  problem  with  the  current  comer 
treatment  derived  from  the  assumption  that  the  entire  tributary  area  of  the  slide  line  which  is 
lumped  at  each  node  effectively  "slides  off'  or  "slides  on1’  to  an  adjacent  node  instantaneously, 
i.e.  the  node  contact  is  binary  (off  or  on)  instead  of  graded  as  the  node’s  tributary  region  slides 
past  the  comer.  This  lumped  behavior  is  fine  tor  many  types  of  problems  but  for  the  jointed 
block  problems  with  the  level  of  discretization  being  used,  the  assumption  of  lumped  contact 
resulted  in  excessive  numerical  vibration. 

In  order  to  correct  for  this  difficulty,  the  computation  of  slide  line  contact  near  coiners 
of  blocks  was  extended  in  order  to  compute  the  proportion  of  a  node’s  tributary  area  which  is 
in  contact  with  the  adjacent  block.  Therefore,  as  a  node  slid  past  the  comer  of  an  adjacent  block 
the  penalty  spring  forces  began  to  decrease  smoothly  in  the  proper  proportion  to  the  remaining 
amount  of  contact  area  instead  of  instantaneously  being  reduced  to  zero.  We  suspect  this 
approach  is  analogous  to  the  usage  by  Itasca  of  "rounded  comers".  While  our  comers  are  not 
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technically  rounded,  the  grading  of  a  penalty  spring’s  resistance  as  nodes  slide  off  of  or  on  to 
a  comer  of  a  block  produces  a  similar  effect. 

The  extensions  to  the  slide  line  logic  were  used  for  the  series  of  larger  benchmark 
problems  helping  to  produce  well-behaved  joint  interaction  computations. 
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